Boring setup stuff
reqpkg <- c("magrittr","dplyr","ggplot2","ggpubr","DT","reshape2")
for (i in reqpkg) {
print(i)
print(packageVersion(i))
suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}
## [1] "magrittr"
## [1] '1.5'
## [1] "dplyr"
## [1] '0.8.5'
## [1] "ggplot2"
## [1] '3.3.0'
## [1] "ggpubr"
## [1] '0.2.5'
## [1] "DT"
## [1] '0.13'
## [1] "reshape2"
## [1] '1.4.3'
reqpkg <- c("purrr","nplr","car","compute.es","effects","multcomp","pastecs","WRS2","outliers")
for (i in reqpkg) {
print(i)
print(packageVersion(i))
suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}
## [1] "purrr"
## [1] '0.3.3'
## [1] "nplr"
## [1] '0.1.7'
## [1] "car"
## [1] '3.0.7'
## [1] "compute.es"
## [1] '0.2.4'
## [1] "effects"
## [1] '4.1.4'
## [1] "multcomp"
## [1] '1.4.12'
## [1] "pastecs"
## [1] '1.3.21'
## [1] "WRS2"
## [1] '1.0.0'
## [1] "outliers"
## [1] '0.14'
theme_set(theme_pubr())
select <- dplyr::select
colors <- c("#FC2A1C", "#103EFB")
uc_colors <- c("#800000", "#D6D6CE")
path <- "../data/R/Obesity/"
# Import data ####
cohort1_1 <- read.csv2(paste0(path, "1_pt1_total.csv"), sep = ",") %>%
set_colnames(paste0("cohort1.1_",colnames(.)))
cohort1_1_dark <- read.csv2(paste0(path, "1_pt1_dark.csv"), sep = ",") %>%
set_colnames(paste0("cohort1.1_",colnames(.)))
cohort1_1_light <- read.csv2(paste0(path, "1_pt1_light.csv"), sep = ",") %>%
set_colnames(paste0("cohort1.1_",colnames(.)))
cohort1_1_covar <- read.csv2(paste0(path, "1_pt1_covar.csv"), sep = ",") %>%
mutate_all(function(x) as.numeric(as.character(x)))
# cohort1_2 <- read.csv2(paste0(path, "1_pt2_total.csv"), sep = ",") %>%
# set_colnames(paste0("cohort1.2_",colnames(.)))
# cohort1_2_dark <- read.csv2(paste0(path, "1_pt2_dark.csv"), sep = ",") %>%
# set_colnames(paste0("cohort1.2_",colnames(.)))
# cohort1_2_light <- read.csv2(paste0(path, "1_pt2_light.csv"), sep = ",") %>%
# set_colnames(paste0("cohort1.2_",colnames(.)))
# cohort1_2_covar <- read.csv2(paste0(path, "1_pt2_covar.csv"), sep = ",") %>%
# mutate_all(function(x) as.numeric(as.character(x)))
cohort1_3 <- read.csv2(paste0(path, "1_pt3_total.csv"), sep = ",") %>%
set_colnames(paste0("cohort1.3_",colnames(.)))
cohort1_3_dark <- read.csv2(paste0(path, "1_pt3_dark.csv"), sep = ",") %>%
set_colnames(paste0("cohort1.3_",colnames(.)))
cohort1_3_light <- read.csv2(paste0(path, "1_pt3_light.csv"), sep = ",") %>%
set_colnames(paste0("cohort1.3_",colnames(.)))
cohort1_3_covar <- read.csv2(paste0(path, "1_pt3_covar.csv"), sep = ",") %>%
mutate_all(function(x) as.numeric(as.character(x)))
cohort2_1 <- read.csv2(paste0(path, "2_pt1_total.csv"), sep = ",") %>%
set_colnames(paste0("cohort2.1_",colnames(.)))
cohort2_1_dark <- read.csv2(paste0(path, "2_pt1_dark.csv"), sep = ",") %>%
set_colnames(paste0("cohort2.1_",colnames(.)))
cohort2_1_light <- read.csv2(paste0(path, "2_pt1_light.csv"), sep = ",") %>%
set_colnames(paste0("cohort2.1_",colnames(.)))
cohort2_1_covar <- read.csv2(paste0(path, "2_pt1_covar.csv"), sep = ",") %>%
mutate_all(function(x) as.numeric(as.character(x)))
cohort2_2 <- read.csv2(paste0(path, "2_pt2_total.csv"), sep = ",") %>%
set_colnames(paste0("cohort2.2_",colnames(.)))
cohort2_2_dark <- read.csv2(paste0(path, "2_pt2_dark.csv"), sep = ",") %>%
set_colnames(paste0("cohort2.2_",colnames(.)))
cohort2_2_light <- read.csv2(paste0(path, "2_pt2_light.csv"), sep = ",") %>%
set_colnames(paste0("cohort2.2_",colnames(.)))
cohort2_2_covar <- read.csv2(paste0(path, "2_pt2_covar.csv"), sep = ",") %>%
mutate_all(function(x) as.numeric(as.character(x)))
cohort2_3 <- read.csv2(paste0(path, "2_pt3_total.csv"), sep = ",") %>%
set_colnames(paste0("cohort2.3_",colnames(.)))
cohort2_3_dark <- read.csv2(paste0(path, "2_pt3_dark.csv"), sep = ",") %>%
set_colnames(paste0("cohort2.3_",colnames(.)))
cohort2_3_light <- read.csv2(paste0(path, "2_pt3_light.csv"), sep = ",") %>%
set_colnames(paste0("cohort2.3_",colnames(.)))
cohort2_3_covar <- read.csv2(paste0(path, "2_pt3_covar.csv"), sep = ",") %>%
mutate_all(function(x) as.numeric(as.character(x)))
cohort3_1 <- read.csv2(paste0(path, "3_pt1_total.csv"), sep = ",") %>%
set_colnames(paste0("cohort3.1_",colnames(.)))
cohort3_1_dark <- read.csv2(paste0(path, "3_pt1_dark.csv"), sep = ",") %>%
set_colnames(paste0("cohort3.1_",colnames(.)))
cohort3_1_light <- read.csv2(paste0(path, "3_pt1_light.csv"), sep = ",") %>%
set_colnames(paste0("cohort3.1_",colnames(.)))
cohort3_1_covar <- read.csv2(paste0(path, "3_pt1_covar.csv"), sep = ",") %>%
mutate_all(function(x) as.numeric(as.character(x)))
cohort3_2 <- read.csv2(paste0(path, "3_pt2_total.csv"), sep = ",") %>%
set_colnames(paste0("cohort3.2_",colnames(.)))
cohort3_2_dark <- read.csv2(paste0(path, "3_pt2_dark.csv"), sep = ",") %>%
set_colnames(paste0("cohort3.2_",colnames(.)))
cohort3_2_light <- read.csv2(paste0(path, "3_pt2_light.csv"), sep = ",") %>%
set_colnames(paste0("cohort3.2_",colnames(.)))
cohort3_2_covar <- read.csv2(paste0(path, "3_pt2_covar.csv"), sep = ",") %>%
mutate_all(function(x) as.numeric(as.character(x)))
cohort3_3 <- read.csv2(paste0(path, "3_pt3_total.csv"), sep = ",") %>%
set_colnames(paste0("cohort3.3_",colnames(.)))
cohort3_3_dark <- read.csv2(paste0(path, "3_pt3_dark.csv"), sep = ",") %>%
set_colnames(paste0("cohort3.3_",colnames(.)))
cohort3_3_light <- read.csv2(paste0(path, "3_pt3_light.csv"), sep = ",") %>%
set_colnames(paste0("cohort3.3_",colnames(.)))
cohort3_3_covar <- read.csv2(paste0(path, "3_pt3_covar.csv"), sep = ",") %>%
mutate_all(function(x) as.numeric(as.character(x)))
cohort2_run1 <- read.csv2(paste0(path, "2_macro13.csv_macro13.csv"), sep = ",") %>%
set_colnames(paste0("cohort2.1_",colnames(.)))
cohort3_run1 <- read.csv2(paste0(path, "3_macro13.csv_macro13.csv"), sep = ",") %>%
set_colnames(paste0("cohort3.1_",colnames(.)))
covar_SPF_pt1 <- c(cohort1_1_covar[[1]][c(2,3)], cohort2_1_covar[[1]][c(6,7,8,9)], cohort3_1_covar[[1]][c(4,5,6)])
covar_GF_pt1 <- c(cohort1_1_covar[[1]][1], cohort2_1_covar[[1]][c(1,2,3,4,5)], cohort3_1_covar[[1]][c(1,2,3)])
covar_SPF_pt3 <- c(cohort1_3_covar[[1]][c(2,3)], cohort2_3_covar[[1]][c(4,5,6,7)], cohort3_3_covar[[1]][c(4,5,6)])
covar_GF_pt3 <- c(cohort1_3_covar[[1]][1], cohort2_3_covar[[1]][c(1,2,3)], cohort3_3_covar[[1]][c(1,2,3)])
total2_pt1 <- cohort1_1 %>%
bind_cols(cohort2_1) %>%
bind_cols(cohort3_1) %>%
select(-contains(c("Enviro", "Time"))) %>%
mutate_all(function(x) as.numeric(as.character(x)))
light_pt1 <- cohort1_1_light %>%
bind_cols(cohort2_1_light) %>%
bind_cols(cohort3_1_light) %>%
select(-contains(c("Enviro", "Time"))) %>%
mutate_all(function(x) as.numeric(as.character(x)))
dark_pt1 <- cohort1_1_dark %>%
bind_cols(cohort2_1_dark) %>%
bind_cols(cohort3_1_dark) %>%
select(-contains(c("Enviro", "Time"))) %>%
mutate_all(function(x) as.numeric(as.character(x)))
total_pt1 <- bind_rows(light_pt1, dark_pt1) %>%
mutate_at(dplyr::vars(contains("VO2")), ~ ifelse(. < 0.6, NA, .))
total2_pt2 <- cohort2_2 %>%
bind_cols(cohort3_2) %>%
select(-contains(c("Enviro", "Time"))) %>%
mutate_all(function(x) as.numeric(as.character(x)))
light_pt2 <- cohort2_2_light %>%
bind_cols(cohort3_2_light) %>%
select(-contains(c("Enviro", "Time"))) %>%
mutate_all(function(x) as.numeric(as.character(x)))
dark_pt2 <- cohort2_2_dark %>%
bind_cols(cohort3_2_light) %>%
select(-contains(c("Enviro", "Time"))) %>%
mutate_all(function(x) as.numeric(as.character(x)))
total_pt2 <- bind_rows(light_pt2, dark_pt2) %>%
mutate_at(dplyr::vars(contains("VO2")), ~ ifelse(. < 0.6, NA, .))
total2_pt3 <- cohort1_3 %>%
bind_cols(cohort2_3) %>%
bind_cols(cohort3_3) %>%
select(-contains(c("Enviro", "Time"))) %>%
mutate_all(function(x) as.numeric(as.character(x)))
light_pt3 <- cohort1_3_light %>%
bind_cols(cohort2_3_light) %>%
bind_cols(cohort3_3_light) %>%
select(-contains(c("Enviro", "Time"))) %>%
mutate_all(function(x) as.numeric(as.character(x)))
dark_pt3 <- cohort1_3_dark %>%
bind_cols(cohort2_3_dark) %>%
bind_cols(cohort3_3_dark) %>%
select(-contains(c("Enviro", "Time"))) %>%
mutate_all(function(x) as.numeric(as.character(x)))
total_pt3 <- bind_rows(light_pt3, dark_pt3)
total_run <- cohort2_run1 %>%
bind_cols(cohort3_run1) %>%
select(-contains(c("Enviro", "Time"))) %>%
mutate_all(function(x) as.numeric(as.character(x)))
# Functions ####
select_GF <- function(df) {
d1 <- df %>% select(contains(c("cohort1.1", "cohort1.2")))
for (name in colnames(d1)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(6,9)) {
d1 %<>% select(-all_of(name))
}
}
d2 <- df %>% select(contains("cohort1.3"))
for (name in colnames(d2)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(11,12)) {
d2 %<>% select(-all_of(name))
}
}
d3 <- df %>% select(contains(c("cohort2.1","cohort2.2")))
for (name in colnames(d3)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(6,9,10,11)) {
d3 %<>% select(-all_of(name))
}
}
d4 <- df %>% select(contains("cohort2.3"))
for (name in colnames(d4)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(9,10,11,12)) {
d4 %<>% select(-all_of(name))
}
}
d5 <- df %>% select(contains(c("cohort3.1","cohort3.2")))
for (name in colnames(d5)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(12,13,14)) {
d5 %<>% select(-all_of(name))
}
}
d6 <- df %>% select(contains("cohort3.3"))
for (name in colnames(d6)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(4,5,6)) {
d6 %<>% select(-all_of(name))
}
}
d <- bind_cols(d1, d2, d3, d4, d5, d6)
return(d)
}
select_SPF <- function(df) {
d1 <- df %>% select(contains(c("cohort1.1","cohort1.2")))
for (name in colnames(d1)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(4,5)) {
d1 %<>% select(-all_of(name))
}
}
d2 <- df %>% select(contains("cohort1.3"))
for (name in colnames(d2)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(10)) {
d2 %<>% select(-all_of(name))
}
}
d3 <- df %>% select(contains(c("cohort2.1","cohort2.2")))
for (name in colnames(d3)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(1,2,3,4,5)) {
d3 %<>% select(-all_of(name))
}
}
d4 <- df %>% select(contains("cohort2.3"))
for (name in colnames(d4)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(1,2,3)) {
d4 %<>% select(-all_of(name))
}
}
d5 <- df %>% select(contains(c("cohort3.1","cohort3.2")))
for (name in colnames(d5)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(9,10,11)) {
d5 %<>% select(-all_of(name))
}
}
d6 <- df %>% select(contains("cohort3.3"))
for (name in colnames(d6)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(1,2,3)) {
d6 %<>% select(-all_of(name))
}
}
d <- bind_cols(d1, d2, d3, d4, d5, d6)
return(d)
}
graph_groups <- function(data, channel, title="", type="m", colors=c("#103EFB", "#FC2A1C")) {
df <- data %>% select(contains(channel))
df.observation <- 1:nrow(df)
SPF <- df %>% select_SPF() %>% set_colnames(paste0("SPF_", colnames(.))) %>% melt(id.vars=NULL) %>% select(-variable)
SPF$variable <- rep("SPF", nrow(SPF))
GF <- df %>% select_GF() %>% set_colnames(paste0("GF_", colnames(.))) %>% melt(id.vars=NULL) %>% select(-variable)
GF$variable <- rep("GF", nrow(GF))
df <- bind_rows(SPF, GF)
df$observation <- rep(df.observation, nrow(df)/length(df.observation))
df <- na.omit(df)
if (type=="m") {
p <- df %>% group_by(variable, observation) %>%
summarise(value=mean(value)) %>%
ggplot(aes(x=observation, y=value, color=variable)) +
geom_rect(aes(xmin=241, xmax=481, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_rect(aes(xmin=721, xmax=961, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_line(size=1) +
scale_color_manual(values = colors) +
scale_x_continuous(breaks=seq(0,nrow(df),240)) +
ggtitle(title) +
ylab(channel)
} else if (type=="s") {
p <- ggplot(df, aes(x=observation, y=value, color=variable)) +
geom_rect(aes(xmin=241, xmax=481, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_rect(aes(xmin=721, xmax=961, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
scale_color_manual(values = colors) +
scale_x_continuous(breaks=seq(0,nrow(df),240)) +
stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", alpha = 0.7) +
ggtitle(title) +
ylab(channel)
}
return(p)
}
graph_mice <- function(data, channel, title="") {
df <- data %>% select(contains(channel))
df.observation <- 1:nrow(df)
SPF <- df %>% select_SPF() %>% set_colnames(paste0("SPF_", colnames(.))) %>% melt()
GF <- df %>% select_GF() %>% set_colnames(paste0("GF_", colnames(.))) %>% melt()
df <- bind_rows(SPF, GF) %>% set_colnames(c("mouse", "channel"))
df$observation <- rep(df.observation, nrow(df)/length(df.observation))
p <- ggplot(df, aes(x=observation, y=channel, color=mouse)) +
geom_rect(aes(xmin=241, xmax=481, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_rect(aes(xmin=721, xmax=961, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
scale_x_continuous(breaks=seq(0,nrow(df),240)) +
geom_line() +
ggtitle(title) +
ylab(channel)
return(p)
}
signif.floor <- function(x, n){
pow <- floor( log10( abs(x) ) ) + 1 - n
y <- floor(x / 10 ^ pow) * 10^pow
# handle the x = 0 case
y[x==0] <- 0
y
}
signif.ceiling <- function(x, n){
pow <- floor( log10( abs(x) ) ) + 1 - n
y <- ceiling(x / 10 ^ pow) * 10^pow
# handle the x = 0 case
y[x==0] <- 0
y
}
find_range <- function(data, channel, segments=10) {
data <- data %>% select(contains(channel)) %>% melt(id.vars = NULL)
r <- range(data$value, na.rm=TRUE)
cut <- max(data$value, na.rm=TRUE)-min(data$value, na.rm=TRUE)
out <- list("range"=r, "cut"=cut/segments)
return(out)
}
rcf <- function(data, channel, id, min, max, b=0.1) {
#' Relative cumulative frequency function
#'
#' Takes a vector and returns a dataframe with relative cumulative frequency per break unit
#'
#' @param data numeric vector
#' @param channel character. Channel to label the rcf output.
#' @param id character. ID for this sample.
#' @param min double. Minimum value of data.
#' @param max double. Maximum value of data.
#' @param b double. Number to cut vector.
#' Before running, check for valid break values to cut the data. One way is to divide range of data by 10 `(max(data)-min(data))/10`.
breaks <- seq(min, max, by=b)
# breaks <- seq(signif.floor(min(data), 2), signif.ceiling(max(data), 2), by=b)
duration.cut <- cut(data, breaks, right = F)
duration.freq <- table(duration.cut)
duration.cumfreq = cumsum(duration.freq) %>% prepend(0)
duration.cumrelfreq = duration.cumfreq/length(data)
out <- duration.cumrelfreq %>% as.data.frame()
out <- cbind(breaks, out) %>% set_colnames(c("breaks", channel)) %>% set_rownames(NULL)
out$id <- rep(id, nrow(out))
return(out)
}
EC50 <- function(x) {
out <- lapply(x, function(y) {
nplr(y$breaks, y[, 2], silent = T, useLog = FALSE)
}) %>%
sapply(function(z) {
getEstimates(z, 0.5)$x
})
out %<>% set_rownames(NULL)
return(out)
}
lm_eqn <- function(df, y, x){
m <- lm(y ~ x, df);
eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2,
list(a = format(unname(coef(m)[1]), digits = 2),
b = format(unname(coef(m)[2]), digits = 2),
r2 = format(summary(m)$r.squared, digits = 3)))
as.character(as.expression(eq));
}
lm_eqn_str <- function(df, y, x){
m <- lm(y ~ x, df);
eq <- sprintf("y = %s + %sx, r^2 = %s", format(unname(coef(m)[1]), digits = 2),format(unname(coef(m)[2]), digits = 2),format(summary(m)$r.squared, digits = 3))
return(eq)
}
ANCOVA <- function(df) {
# Takes dataframe with column 1 as independent variable (e.g., WT), column 2 as dependent variable, and column 3 as covariate
cat("linear fit of dependent variable with covariates\n")
p1 <- ggplot(df, aes_string(x=names(df)[2], y=names(df)[3])) +
geom_smooth(method = "lm", se = FALSE) +
geom_point() +
# geom_text(x = mean(na.omit(df[[2]])), y = mean(na.omit(df[[3]])), label = lm_eqn(df, df[[3]], df[[2]]), parse = TRUE) +
geom_text(x = mean(na.omit(df[[2]])), y = (mean(na.omit(df[[3]])) + IQR(df[[3]], na.rm = T)), label = lm_eqn(df, df[[3]], df[[2]]), parse = TRUE) +
ggtitle("Linear model - DV and Covar") +
xlab("DV") + ylab("Covar")
print(p1)
m <- lm(df[[3]] ~ df[[2]], df)
summary(m) %>% print()
lm_eqn_str(df, df[[3]], df[[2]]) %>% cat()
plot1 <- ggplot(df, aes_string(x=names(df)[1], y=names(df)[2])) +
geom_boxplot() +
# ggtitle("DV vs IV") +
xlab("IV") + ylab("DV")
plot2 <- ggplot(df, aes_string(x=names(df)[1], y=names(df)[3])) +
geom_boxplot() +
# ggtitle("Covar vs IV") +
xlab("IV") + ylab("Covar")
p2 <- ggarrange(plot1, plot2,
labels = c("DV vs IV", "Covar vs IV"),
ncol = 2)
# grid.arrange(plot1, plot2, ncol=2)
print(p2)
cat("\n―――――――――――――――――――――\n")
cat("Some stats on the variables\n")
cat("DV - IV:\n")
by(df[[2]], df[[1]], stat.desc) %>% print()
cat("Covar - IV:\n")
by(df[[3]], df[[1]], stat.desc) %>% print()
cat("levene's test\n")
leveneTest(df[[2]], df[[1]], center = median) %>% print()
cat("―――――――――――――――――――――\n")
cat("ANOVA of DV with IV\n")
model <- aov(df[[2]] ~ df[[1]], df)
summary(model) %>% print()
cat("―――――――――――――――――――――\n")
cat("ANOVA of covariate with IV, independent of DV\n")
model <- aov(df[[3]] ~ df[[1]], df)
summary(model) %>% print()
cat("―――――――――――――――――――――\n")
cat("ANOVA of DV with IV, normalized to covariate\n")
c <- df[[2]]/df[[3]]
model <- aov(c ~ df[[1]], df)
summary(model) %>% print()
cat("―――――――――――――――――――――\n")
cat("ANCOVA\n")
model <- aov(df[[2]] ~ df[[1]] + df[[3]], data = df)
summary(model) %>% print()
out <- Anova(model, type= "III")
print(out)
return(out)
}
check_ANCOVA <- function(res) {
if (res$`Pr(>F)`[2] < 0.05 && res$`Pr(>F)`[3] < 0.05) {
print(paste("Varies significantly with DV and Covar,"), res$`Pr(>F)`[2], res$`Pr(>F)`[3], sep=" ")
} else if (res$`Pr(>F)`[2] < 0.05 && res$`Pr(>F)`[3] > 0.05) {
print(paste("Varies significantly with DV,", "p =", res$`Pr(>F)`[2], sep =" "))
} else if (res$`Pr(>F)`[2] > 0.05 && res$`Pr(>F)`[3] < 0.05) {
print(paste("Varies significantly with Covar,", "p =", res$`Pr(>F)`[3], sep=" "))
}
}
rcf_pipe <- function(data, cycle, channel, segments = 50, id1 = "SPF", id2 = "GF", colors = c("#128D15", "#FD9226")) {
r <- find_range(data, channel, segments)
print("Range:")
print(r)
data <- data %>% select(contains(channel))
d <- list()
d[[id1]] <- data %>% select_SPF()
d[[id2]] <- data %>% select_GF()
rcf_res <- list()
rcf_res[[id1]] <- lapply(d[[id1]], function(x) rcf(x, channel, id1, r$range[1], r$range[2], r$cut))
rcf_res[[id2]] <- lapply(d[[id2]], function(x) rcf(x, channel, id2, r$range[1], r$range[2], r$cut))
data <- bind_rows(bind_rows(rcf_res[[id1]]), bind_rows(rcf_res[[id2]])) %>% split(.$id)
Models <- lapply(data, function(x){
nplr(x$breaks, x[[channel]], silent = TRUE, useLog = F)
})
overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main=paste0("RCF of ", id1, " vs ", id2, " mice ", channel, ", ", cycle), cex.main=1.5, lwd = 3, Cols = c("#103EFB", "#FC2A1C"))
abline(h=0.5)
print(paste0(id1, " EC50: ", getEstimates(Models[[id1]], 0.5)$x))
print(paste0(id2, " EC50: ", getEstimates(Models[[id2]], 0.5)$x))
return(list(rcf_res[[id1]], rcf_res[[id2]]))
}
VO2
# VO2 ####
cuts <- seq(0,nrow(total_run),240)
cuts <- sapply(cuts, function(x) x+1)
channel = "VO2"
p <- graph_groups(total_run, channel, title = paste0(channel, " before and after adding HF diet"))
p + geom_rect(aes(xmin=cuts[6], xmax=cuts[7], ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_rect(aes(xmin=cuts[8], xmax=cuts[9], ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_rect(aes(xmin=cuts[10], xmax=cuts[11], ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_rect(aes(xmin=cuts[12], xmax=cuts[13], ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_rect(aes(xmin=cuts[14], xmax=cuts[15], ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_line(size=1) +
geom_vline(xintercept = 1500)
# p1 <- graph_groups(total2_pt1, channel, title = paste0(channel, " before HF diet"))
# p2 <- graph_groups(total2_pt2, channel, title = paste0(channel, " immediately after HF diet"))
# p3 <- graph_groups(total2_pt3, channel, title = paste0(channel, " 1 month after HF diet"))
p1 <- graph_groups(total2_pt1, channel)
p2 <- graph_groups(total2_pt2, channel)
p3 <- graph_groups(total2_pt3, channel)
leg <- get_legend(p1)
p1 <- p1 + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
p2 <- p2 + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
p3 <- p3 + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
f <- ggarrange(p1, p2, p3, leg)
annotate_figure(f, top = text_grob(paste0(channel, ", before, after, and a month after adding HF diet"), face = "bold"), left = text_grob(channel), bottom = text_grob("observation"))
Before HF diet
# pt1 VO2 ####
pt = "before HF diet"
cycle = "Total"
total_pt1.VO2.rcf <- rcf_pipe(total_pt1, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.8746158 4.9918380
##
## $cut
## [1] 0.08234444
## [1] "SPF EC50: 1.75001617929315"
## [1] "GF EC50: 1.60877208757266"
t <- total_pt1.VO2.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)
ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50", pt, sep = " "))
dixon.test(t[[1]])
##
## Dixon test for outliers
##
## data: t[[1]]
## Q.cohort2.1_VO2_M_6 = 0.42252, p-value = 0.2344
## alternative hypothesis: highest value 2.08865368747895 is an outlier
dixon.test(t[[1]][-which.max(t[[1]])])
##
## Dixon test for outliers
##
## data: t[[1]][-which.max(t[[1]])]
## Q.cohort1.1_VO2_M_9 = 0.43748, p-value = 0.2774
## alternative hypothesis: lowest value 1.46653571141594 is an outlier
dixon.test(t[[2]])
##
## Dixon test for outliers
##
## data: t[[2]]
## Q.cohort3.1_VO2_M_10 = 0.049739, p-value = 0.2924
## alternative hypothesis: lowest value 1.42989814868382 is an outlier
f <- ggarrange(
ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", sep = " "), face = "bold"))
shapiro.test(t[[1]])
##
## Shapiro-Wilk normality test
##
## data: t[[1]]
## W = 0.97779, p-value = 0.952
shapiro.test(t[[2]])
##
## Shapiro-Wilk normality test
##
## data: t[[2]]
## W = 0.91597, p-value = 0.3599
var.test(t[[1]], t[[2]])
##
## F test to compare two variances
##
## data: t[[1]] and t[[2]]
## F = 2.3404, num df = 8, denom df = 8, p-value = 0.2504
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.5279073 10.3753815
## sample estimates:
## ratio of variances
## 2.34035
# t.test(t[[1]], t[[2]])
# If var.test is not significant, set var.equal to true!
t.test(t[[1]], t[[2]], var.equal = T)
##
## Two Sample t-test
##
## data: t[[1]] and t[[2]]
## t = 2.0526, df = 16, p-value = 0.05684
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.00474476 0.29423557
## sample estimates:
## mean of x mean of y
## 1.749973 1.605228
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.7152 -1.5104 -0.0060 0.7461 5.6158
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 27.9665 5.7561 4.859 0.000174 ***
## df[[2]] 0.8346 3.4159 0.244 0.810083
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.297 on 16 degrees of freedom
## Multiple R-squared: 0.003717, Adjusted R-squared: -0.05855
## F-statistic: 0.0597 on 1 and 16 DF, p-value: 0.8101
##
## y = 28 + 0.83x, r^2 = 0.00372
##
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 1.42989815 1.75896428 0.32906613
## sum median mean SE.mean CI.mean.0.95 var
## 14.44704950 1.63768933 1.60522772 0.03858337 0.08897341 0.01339809
## std.dev coef.var
## 0.11575011 0.07210822
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 1.46653571 2.08865369 0.62211798
## sum median mean SE.mean CI.mean.0.95 var
## 15.74975817 1.75334199 1.74997313 0.05902563 0.13611334 0.03135622
## std.dev coef.var
## 0.17707688 0.10118834
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 29.00000000 35.00000000 6.00000000
## sum median mean SE.mean CI.mean.0.95 var
## 276.00000000 30.00000000 30.66666667 0.68211273 1.57295478 4.18750000
## std.dev coef.var
## 2.04633819 0.06672842
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 25.60000000 30.40000000 4.80000000
## sum median mean SE.mean CI.mean.0.95 var
## 252.60000000 27.80000000 28.06666667 0.53800041 1.24063118 2.60500000
## std.dev coef.var
## 1.61400124 0.05750598
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.8135 0.3805
## 16
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0943 0.09428 4.213 0.0568 .
## Residuals 16 0.3580 0.02238
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 30.42 30.420 8.957 0.00861 **
## Residuals 16 54.34 3.396
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0004462 0.0004462 20.34 0.000356 ***
## Residuals 16 0.0003509 0.0000219
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.09428 0.09428 5.067 0.0398 *
## df[[3]] 1 0.07893 0.07893 4.242 0.0572 .
## Residuals 15 0.27910 0.01861
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.010936 1 0.5878 0.455182
## df[[1]] 0.171531 1 9.2187 0.008334 **
## df[[3]] 0.078931 1 4.2421 0.057229 .
## Residuals 0.279103 15
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.00833428507313233"
cycle = "Light"
light_pt1.VO2.rcf <- rcf_pipe(light_pt1, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.8746158 2.9373260
##
## $cut
## [1] 0.0412542
## [1] "SPF EC50: 1.58906026308396"
## [1] "GF EC50: 1.44574161742324"
t <- light_pt1.VO2.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)
ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50", pt, sep = " "))
dixon.test(t[[1]])
##
## Dixon test for outliers
##
## data: t[[1]]
## Q.cohort2.1_VO2_M_6 = 0.19725, p-value = 0.9629
## alternative hypothesis: highest value 1.79542307671054 is an outlier
dixon.test(t[[2]])
##
## Dixon test for outliers
##
## data: t[[2]]
## Q.cohort1.1_VO2_M_4 = 0.076364, p-value = 0.4349
## alternative hypothesis: lowest value 1.26352543847839 is an outlier
f <- ggarrange(
ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))
shapiro.test(t[[1]])
##
## Shapiro-Wilk normality test
##
## data: t[[1]]
## W = 0.92042, p-value = 0.3958
shapiro.test(t[[2]])
##
## Shapiro-Wilk normality test
##
## data: t[[2]]
## W = 0.82751, p-value = 0.04187
var.test(t[[1]], t[[2]])
##
## F test to compare two variances
##
## data: t[[1]] and t[[2]]
## F = 2.499, num df = 8, denom df = 8, p-value = 0.2167
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.5636947 11.0787387
## sample estimates:
## ratio of variances
## 2.499005
t.test(t[[1]], t[[2]], var.equal = T)
##
## Two Sample t-test
##
## data: t[[1]] and t[[2]]
## t = 2.2009, df = 16, p-value = 0.04277
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.005185257 0.276681933
## sample estimates:
## mean of x mean of y
## 1.568389 1.427456
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.7341 -1.5320 -0.0815 0.8395 5.6255
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.913 5.584 5.177 9.17e-05 ***
## df[[2]] 0.303 3.711 0.082 0.936
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.301 on 16 degrees of freedom
## Multiple R-squared: 0.0004167, Adjusted R-squared: -0.06206
## F-statistic: 0.00667 on 1 and 16 DF, p-value: 0.9359
##
## y = 29 + 0.3x, r^2 = 0.000417
##
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 1.26352544 1.52389696 0.26037152
## sum median mean SE.mean CI.mean.0.95 var
## 12.84710286 1.48944874 1.42745587 0.03423306 0.07894159 0.01054712
## std.dev coef.var
## 0.10269919 0.07194562
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 1.34191921 1.79542308 0.45350387
## sum median mean SE.mean CI.mean.0.95 var
## 14.11550522 1.63088754 1.56838947 0.05411646 0.12479277 0.02635732
## std.dev coef.var
## 0.16234937 0.10351343
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 29.00000000 35.00000000 6.00000000
## sum median mean SE.mean CI.mean.0.95 var
## 276.00000000 30.00000000 30.66666667 0.68211273 1.57295478 4.18750000
## std.dev coef.var
## 2.04633819 0.06672842
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 25.60000000 30.40000000 4.80000000
## sum median mean SE.mean CI.mean.0.95 var
## 252.60000000 27.80000000 28.06666667 0.53800041 1.24063118 2.60500000
## std.dev coef.var
## 1.61400124 0.05750598
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.482 0.2411
## 16
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.08938 0.08938 4.844 0.0428 *
## Residuals 16 0.29524 0.01845
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 30.42 30.420 8.957 0.00861 **
## Residuals 16 54.34 3.396
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0003909 0.0003909 19.97 0.000387 ***
## Residuals 16 0.0003131 0.0000196
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.08938 0.08938 5.636 0.0314 *
## df[[3]] 1 0.05736 0.05736 3.617 0.0766 .
## Residuals 15 0.23788 0.01586
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.010671 1 0.6729 0.424895
## df[[1]] 0.146579 1 9.2430 0.008267 **
## df[[3]] 0.057359 1 3.6169 0.076576 .
## Residuals 0.237876 15
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.00826660982644304"
cycle = "Dark"
dark_pt1.VO2.rcf <- rcf_pipe(dark_pt1, cycle, channel)
## [1] "Range:"
## $range
## [1] 1.000101 4.991838
##
## $cut
## [1] 0.07983474
## [1] "SPF EC50: 1.9341144311299"
## [1] "GF EC50: 1.83403636087834"
t <- dark_pt1.VO2.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)
ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50", pt, sep = " "))
dixon.test(t[[1]])
##
## Dixon test for outliers
##
## data: t[[1]]
## Q.cohort2.1_VO2_M_6 = 0.47365, p-value = 0.1482
## alternative hypothesis: highest value 2.37415711632301 is an outlier
dixon.test(t[[2]])
##
## Dixon test for outliers
##
## data: t[[2]]
## Q.cohort3.1_VO2_M_10 = 0.097525, p-value = 0.5498
## alternative hypothesis: lowest value 1.59851055132634 is an outlier
f <- ggarrange(
ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))
shapiro.test(t[[1]])
##
## Shapiro-Wilk normality test
##
## data: t[[1]]
## W = 0.94522, p-value = 0.6378
shapiro.test(t[[2]])
##
## Shapiro-Wilk normality test
##
## data: t[[2]]
## W = 0.93588, p-value = 0.5393
var.test(t[[1]], t[[2]])
##
## F test to compare two variances
##
## data: t[[1]] and t[[2]]
## F = 2.1262, num df = 8, denom df = 8, p-value = 0.3065
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.4796022 9.4260035
## sample estimates:
## ratio of variances
## 2.126201
t.test(t[[1]], t[[2]], var.equal = T)
##
## Two Sample t-test
##
## data: t[[1]] and t[[2]]
## t = 1.6297, df = 16, p-value = 0.1227
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.04176041 0.31941330
## sample estimates:
## mean of x mean of y
## 1.953008 1.814181
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6744 -1.3941 0.0982 0.7798 5.6352
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.732 5.541 4.824 0.000187 ***
## df[[2]] 1.399 2.928 0.478 0.639377
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.285 on 16 degrees of freedom
## Multiple R-squared: 0.01406, Adjusted R-squared: -0.04756
## F-statistic: 0.2281 on 1 and 16 DF, p-value: 0.6394
##
## y = 27 + 1.4x, r^2 = 0.0141
##
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 1.59851055 2.00713925 0.40862869
## sum median mean SE.mean CI.mean.0.95 var
## 16.32762964 1.82944446 1.81418107 0.04817938 0.11110185 0.02089127
## std.dev coef.var
## 0.14453814 0.07967129
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 1.60328120 2.37415712 0.77087591
## sum median mean SE.mean CI.mean.0.95 var
## 17.57706766 1.93427500 1.95300752 0.07025276 0.16200316 0.04441905
## std.dev coef.var
## 0.21075829 0.10791473
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 29.00000000 35.00000000 6.00000000
## sum median mean SE.mean CI.mean.0.95 var
## 276.00000000 30.00000000 30.66666667 0.68211273 1.57295478 4.18750000
## std.dev coef.var
## 2.04633819 0.06672842
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 25.60000000 30.40000000 4.80000000
## sum median mean SE.mean CI.mean.0.95 var
## 252.60000000 27.80000000 28.06666667 0.53800041 1.24063118 2.60500000
## std.dev coef.var
## 1.61400124 0.05750598
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.2463 0.6264
## 16
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0867 0.08673 2.656 0.123
## Residuals 16 0.5225 0.03266
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 30.42 30.420 8.957 0.00861 **
## Residuals 16 54.34 3.396
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0004835 0.0004835 15.65 0.00113 **
## Residuals 16 0.0004944 0.0000309
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0867 0.08673 3.176 0.0950 .
## df[[3]] 1 0.1128 0.11284 4.132 0.0602 .
## Residuals 15 0.4096 0.02731
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.00997 1 0.3651 0.55474
## df[[1]] 0.19101 1 6.9941 0.01839 *
## df[[3]] 0.11284 1 4.1320 0.06017 .
## Residuals 0.40964 15
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.0183871212047179"
Immediately after HF diet
# pt2 VO2 ####
# delete these once I add cohort1_2!!
covar_SPF_pt1 <- c(cohort2_1_covar[[1]][c(6,7,8,9)], cohort3_1_covar[[1]][c(4,5,6)])
covar_GF_pt1 <- c(cohort2_1_covar[[1]][c(1,2,3,4,5)], cohort3_1_covar[[1]][c(1,2,3)])
pt = "immediately after HF diet"
cycle = "Total"
total_pt2.VO2.rcf <- rcf_pipe(total_pt2, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.8753851 3.3065780
##
## $cut
## [1] 0.04862386
## [1] "SPF EC50: 1.81300642705308"
## [1] "GF EC50: 1.6968627087087"
t <- total_pt2.VO2.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)
ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50", pt, sep = " "))
dixon.test(t[[1]])
##
## Dixon test for outliers
##
## data: t[[1]]
## Q.cohort2.2_VO2_M_6 = 0.39782, p-value = 0.2693
## alternative hypothesis: highest value 2.33215178662138 is an outlier
dixon.test(t[[2]])
##
## Dixon test for outliers
##
## data: t[[2]]
## Q.cohort2.2_VO2_M_1 = 0.39801, p-value = 0.3661
## alternative hypothesis: lowest value 1.46610566358194 is an outlier
f <- ggarrange(
ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))
shapiro.test(t[[1]])
##
## Shapiro-Wilk normality test
##
## data: t[[1]]
## W = 0.94242, p-value = 0.6606
shapiro.test(t[[2]])
##
## Shapiro-Wilk normality test
##
## data: t[[2]]
## W = 0.94787, p-value = 0.6897
var.test(t[[1]], t[[2]])
##
## F test to compare two variances
##
## data: t[[1]] and t[[2]]
## F = 5.7567, num df = 6, denom df = 7, p-value = 0.03675
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 1.124668 32.787226
## sample estimates:
## ratio of variances
## 5.756719
t.test(t[[1]], t[[2]], var.equal = FALSE) # F-test p-val: 0.03675
##
## Welch Two Sample t-test
##
## data: t[[1]] and t[[2]]
## t = 1.2272, df = 7.808, p-value = 0.2555
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.1244126 0.4048969
## sample estimates:
## mean of x mean of y
## 1.829071 1.688829
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6643 -1.2707 -0.2855 0.5268 5.1774
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 27.725 4.901 5.657 7.84e-05 ***
## df[[2]] 1.141 2.774 0.411 0.688
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.23 on 13 degrees of freedom
## Multiple R-squared: 0.01285, Adjusted R-squared: -0.06309
## F-statistic: 0.1692 on 1 and 13 DF, p-value: 0.6875
##
## y = 28 + 1.1x, r^2 = 0.0128
##
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 8.00000000 0.00000000 0.00000000 1.46610566 1.83829826 0.37219260
## sum median mean SE.mean CI.mean.0.95 var
## 13.51063076 1.70811396 1.68882884 0.04150956 0.09815452 0.01378435
## std.dev coef.var
## 0.11740677 0.06951964
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 7.00000000 0.00000000 0.00000000 1.40283486 2.33215179 0.92931693
## sum median mean SE.mean CI.mean.0.95 var
## 12.80349695 1.79453176 1.82907099 0.10647108 0.26052535 0.07935264
## std.dev coef.var
## 0.28169600 0.15401043
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 8.00000000 0.00000000 0.00000000 29.00000000 35.00000000 6.00000000
## sum median mean SE.mean CI.mean.0.95 var
## 246.50000000 30.00000000 30.81250000 0.75555975 1.78661491 4.56696429
## std.dev coef.var
## 2.13704569 0.06935645
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 7.0000000 0.0000000 0.0000000 26.3000000 30.4000000 4.1000000
## sum median mean SE.mean CI.mean.0.95 var
## 199.4000000 28.1000000 28.4857143 0.5629127 1.3773978 2.2180952
## std.dev coef.var
## 1.4893271 0.0522833
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.8247 0.1998
## 13
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0734 0.07343 1.667 0.219
## Residuals 13 0.5726 0.04405
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 20.21 20.212 5.803 0.0315 *
## Residuals 13 45.28 3.483
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0003232 0.0003232 7.064 0.0197 *
## Residuals 13 0.0005949 0.0000458
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0734 0.07343 1.805 0.204
## df[[3]] 1 0.0845 0.08446 2.076 0.175
## Residuals 12 0.4881 0.04068
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.00608 1 0.1494 0.70590
## df[[1]] 0.14959 1 3.6773 0.07927 .
## df[[3]] 0.08446 1 2.0763 0.17518
## Residuals 0.48814 12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
cycle = "Light"
light_pt2.VO2.rcf <- rcf_pipe(light_pt2, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.1013106 3.2033890
##
## $cut
## [1] 0.06204157
## [1] "SPF EC50: 1.71309371902381"
## [1] "GF EC50: 1.59241129302527"
t <- light_pt2.VO2.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)
ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50", pt, sep = " "))
dixon.test(t[[1]])
##
## Dixon test for outliers
##
## data: t[[1]]
## Q.cohort2.2_VO2_M_6 = 0.40248, p-value = 0.2596
## alternative hypothesis: highest value 2.05964348633281 is an outlier
dixon.test(t[[2]])
##
## Dixon test for outliers
##
## data: t[[2]]
## Q.cohort3.2_VO2_M_10 = 0.25228, p-value = 0.8294
## alternative hypothesis: lowest value 1.39075334998537 is an outlier
f <- ggarrange(
ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))
shapiro.test(t[[1]])
##
## Shapiro-Wilk normality test
##
## data: t[[1]]
## W = 0.95817, p-value = 0.8029
shapiro.test(t[[2]])
##
## Shapiro-Wilk normality test
##
## data: t[[2]]
## W = 0.91587, p-value = 0.3973
var.test(t[[1]], t[[2]])
##
## F test to compare two variances
##
## data: t[[1]] and t[[2]]
## F = 2.8859, num df = 6, denom df = 7, p-value = 0.1917
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.5638061 16.4365354
## sample estimates:
## ratio of variances
## 2.885896
t.test(t[[1]], t[[2]], var.equal = T)
##
## Two Sample t-test
##
## data: t[[1]] and t[[2]]
## t = 1.5625, df = 13, p-value = 0.1422
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.04971928 0.30961818
## sample estimates:
## mean of x mean of y
## 1.713522 1.583573
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5318 -1.2153 -0.2414 0.5103 5.2058
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.068 5.855 4.794 0.000351 ***
## df[[2]] 1.009 3.543 0.285 0.780371
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.238 on 13 degrees of freedom
## Multiple R-squared: 0.006196, Adjusted R-squared: -0.07025
## F-statistic: 0.08104 on 1 and 13 DF, p-value: 0.7804
##
## y = 28 + 1x, r^2 = 0.0062
##
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 8.00000000 0.00000000 0.00000000 1.39075335 1.71114273 0.32038938
## sum median mean SE.mean CI.mean.0.95 var
## 12.66858433 1.60455785 1.58357304 0.04154122 0.09822937 0.01380538
## std.dev coef.var
## 0.11749630 0.07419696
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 7.00000000 0.00000000 0.00000000 1.40092810 2.05964349 0.65871539
## sum median mean SE.mean CI.mean.0.95 var
## 11.99465742 1.72520434 1.71352249 0.07544240 0.18460091 0.03984089
## std.dev coef.var
## 0.19960184 0.11648627
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 8.00000000 0.00000000 0.00000000 29.00000000 35.00000000 6.00000000
## sum median mean SE.mean CI.mean.0.95 var
## 246.50000000 30.00000000 30.81250000 0.75555975 1.78661491 4.56696429
## std.dev coef.var
## 2.13704569 0.06935645
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 7.0000000 0.0000000 0.0000000 26.3000000 30.4000000 4.1000000
## sum median mean SE.mean CI.mean.0.95 var
## 199.4000000 28.1000000 28.4857143 0.5629127 1.3773978 2.2180952
## std.dev coef.var
## 1.4893271 0.0522833
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.4183 0.529
## 13
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0630 0.06304 2.442 0.142
## Residuals 13 0.3357 0.02582
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 20.21 20.212 5.803 0.0315 *
## Residuals 13 45.28 3.483
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0002788 2.788e-04 10.46 0.00653 **
## Residuals 13 0.0003466 2.666e-05
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.06304 0.06304 2.665 0.129
## df[[3]] 1 0.05177 0.05177 2.188 0.165
## Residuals 12 0.28391 0.02366
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.013909 1 0.5879 0.45806
## df[[1]] 0.112346 1 4.7485 0.04997 *
## df[[3]] 0.051772 1 2.1882 0.16483
## Residuals 0.283911 12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.0499736001563136"
cycle = "Dark"
dark_pt2.VO2.rcf <- rcf_pipe(dark_pt2, cycle, channel)
## [1] "Range:"
## $range
## [1] 1.032885 3.306578
##
## $cut
## [1] 0.04547386
## [1] "SPF EC50: 1.92194220183729"
## [1] "GF EC50: 1.84062090652473"
t <- dark_pt2.VO2.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)
ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50", pt, sep = " "))
dixon.test(t[[1]])
##
## Dixon test for outliers
##
## data: t[[1]]
## Q.cohort2.2_VO2_M_6 = 0.31187, p-value = 0.4956
## alternative hypothesis: highest value 2.60800348939703 is an outlier
dixon.test(t[[2]])
##
## Dixon test for outliers
##
## data: t[[2]]
## Q.cohort3.2_VO2_M_9 = 0.37516, p-value = 0.4268
## alternative hypothesis: lowest value 1.46852607705669 is an outlier
f <- ggarrange(
ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))
shapiro.test(t[[1]])
##
## Shapiro-Wilk normality test
##
## data: t[[1]]
## W = 0.97884, p-value = 0.9537
shapiro.test(t[[2]])
##
## Shapiro-Wilk normality test
##
## data: t[[2]]
## W = 0.84905, p-value = 0.09319
var.test(t[[1]], t[[2]])
##
## F test to compare two variances
##
## data: t[[1]] and t[[2]]
## F = 4.5986, num df = 6, denom df = 7, p-value = 0.06562
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.898416 26.191357
## sample estimates:
## ratio of variances
## 4.598629
t.test(t[[1]], t[[2]], var.equal = T)
##
## Two Sample t-test
##
## data: t[[1]] and t[[2]]
## t = 0.93562, df = 13, p-value = 0.3665
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.1880892 0.4754632
## sample estimates:
## mean of x mean of y
## 1.955075 1.811388
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6876 -1.2898 -0.2909 0.5627 5.2043
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.3373 3.8384 7.383 5.33e-06 ***
## df[[2]] 0.7396 2.0202 0.366 0.72
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.233 on 13 degrees of freedom
## Multiple R-squared: 0.01021, Adjusted R-squared: -0.06593
## F-statistic: 0.134 on 1 and 13 DF, p-value: 0.7202
##
## y = 28 + 0.74x, r^2 = 0.0102
##
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 8.00000000 0.00000000 0.00000000 1.46852608 1.97171082 0.50318475
## sum median mean SE.mean CI.mean.0.95 var
## 14.49110787 1.89549878 1.81138848 0.06431407 0.15207860 0.03309039
## std.dev coef.var
## 0.18190765 0.10042443
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 7.0000000 0.0000000 0.0000000 1.3997469 2.6080035 1.2082566
## sum median mean SE.mean CI.mean.0.95 var
## 13.6855283 1.8564543 1.9550755 0.1474403 0.3607734 0.1521705
## std.dev coef.var
## 0.3900903 0.1995270
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 8.00000000 0.00000000 0.00000000 29.00000000 35.00000000 6.00000000
## sum median mean SE.mean CI.mean.0.95 var
## 246.50000000 30.00000000 30.81250000 0.75555975 1.78661491 4.56696429
## std.dev coef.var
## 2.13704569 0.06935645
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 7.0000000 0.0000000 0.0000000 26.3000000 30.4000000 4.1000000
## sum median mean SE.mean CI.mean.0.95 var
## 199.4000000 28.1000000 28.4857143 0.5629127 1.3773978 2.2180952
## std.dev coef.var
## 1.4893271 0.0522833
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.9939 0.1814
## 13
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0771 0.07708 0.875 0.367
## Residuals 13 1.1447 0.08805
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 20.21 20.212 5.803 0.0315 *
## Residuals 13 45.28 3.483
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0003571 0.0003571 3.648 0.0785 .
## Residuals 13 0.0012727 0.0000979
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0771 0.07708 0.887 0.365
## df[[3]] 1 0.1023 0.10226 1.177 0.299
## Residuals 12 1.0424 0.08687
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.00571 1 0.0657 0.8020
## df[[1]] 0.16687 1 1.9210 0.1910
## df[[3]] 0.10226 1 1.1773 0.2992
## Residuals 1.04239 12
check_ANCOVA(res)
1 month after HF diet
# pt3 VO2 ####
# cohort2.3_VO2_M_1 is an outlier.
total_pt3 %<>% select(-cohort2.3_VO2_M_1)
light_pt3 %<>% select(-cohort2.3_VO2_M_1)
dark_pt3 %<>% select(-cohort2.3_VO2_M_1)
covar_SPF_pt3 <- c(cohort1_3_covar[[1]][c(2,3)], cohort2_3_covar[[1]][c(4,5,6,7)], cohort3_3_covar[[1]][c(4,5,6)])
covar_GF_pt3 <- c(cohort1_3_covar[[1]][1], cohort2_3_covar[[1]][c(2,3)], cohort3_3_covar[[1]][c(1,2,3)])
pt = "1 month after HF diet"
cycle = "Total"
total_pt3.VO2.rcf <- rcf_pipe(total_pt3, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.1282987 3.6305530
##
## $cut
## [1] 0.07004509
## [1] "SPF EC50: 1.86394965004008"
## [1] "GF EC50: 1.68055953099069"
t <- total_pt3.VO2.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt3,covar_GF_pt3)) %>% set_rownames(NULL)
ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50", pt, sep = " "))
dixon.test(t[[1]])
##
## Dixon test for outliers
##
## data: t[[1]]
## Q.cohort2.3_VO2_M_9 = 0.13836, p-value = 0.7538
## alternative hypothesis: lowest value 1.64628551511344 is an outlier
dixon.test(t[[2]])
##
## Dixon test for outliers
##
## data: t[[2]]
## Q.cohort3.3_VO2_M_3 = 0.31403, p-value = 0.6131
## alternative hypothesis: highest value 1.7862408128281 is an outlier
f <- ggarrange(
ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))
shapiro.test(t[[1]])
##
## Shapiro-Wilk normality test
##
## data: t[[1]]
## W = 0.95564, p-value = 0.7518
shapiro.test(t[[2]])
##
## Shapiro-Wilk normality test
##
## data: t[[2]]
## W = 0.99256, p-value = 0.9945
var.test(t[[1]], t[[2]])
##
## F test to compare two variances
##
## data: t[[1]] and t[[2]]
## F = 6.0101, num df = 8, denom df = 5, p-value = 0.06389
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.889446 28.952498
## sample estimates:
## ratio of variances
## 6.010139
t.test(t[[1]], t[[2]], var.equal = T)
##
## Two Sample t-test
##
## data: t[[1]] and t[[2]]
## t = 2.758, df = 13, p-value = 0.01629
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.0401974 0.3308189
## sample estimates:
## mean of x mean of y
## 1.875537 1.690029
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.7604 -2.7232 -0.9105 2.5311 5.9239
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.707 9.858 3.115 0.00821 **
## df[[2]] 1.014 5.454 0.186 0.85538
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.16 on 13 degrees of freedom
## Multiple R-squared: 0.002652, Adjusted R-squared: -0.07407
## F-statistic: 0.03457 on 1 and 13 DF, p-value: 0.8554
##
## y = 31 + 1x, r^2 = 0.00265
##
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 6.000000000 0.000000000 0.000000000 1.605187543 1.786240813 0.181053270
## sum median mean SE.mean CI.mean.0.95 var
## 10.140175624 1.683826709 1.690029271 0.025783815 0.066279407 0.003988831
## std.dev coef.var
## 0.063157191 0.037370472
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 1.64628552 2.09318823 0.44690271
## sum median mean SE.mean CI.mean.0.95 var
## 16.87983684 1.87746493 1.87553743 0.05161118 0.11901560 0.02397343
## std.dev coef.var
## 0.15483355 0.08255423
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 6.00000000 0.00000000 0.00000000 28.70000000 33.20000000 4.50000000
## sum median mean SE.mean CI.mean.0.95 var
## 179.60000000 29.35000000 29.93333333 0.67806915 1.74303225 2.75866667
## std.dev coef.var
## 1.66092344 0.05548742
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 30.90000000 38.30000000 7.40000000
## sum median mean SE.mean CI.mean.0.95 var
## 308.40000000 34.90000000 34.26666667 0.82276634 1.89730257 6.09250000
## std.dev coef.var
## 2.46829901 0.07203207
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 3.6773 0.07739 .
## 13
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.1239 0.12389 7.607 0.0163 *
## Residuals 13 0.2117 0.01629
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 67.60 67.60 14.05 0.00243 **
## Residuals 13 62.53 4.81
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0000070 7.030e-06 0.151 0.704
## Residuals 13 0.0006059 4.661e-05
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.1239 0.12389 13.84 0.00293 **
## df[[3]] 1 0.1043 0.10428 11.65 0.00515 **
## Residuals 12 0.1075 0.00895
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.58516 1 65.349 3.379e-06 ***
## df[[1]] 0.22728 1 25.382 0.0002903 ***
## df[[3]] 0.10428 1 11.646 0.0051489 **
## Residuals 0.10745 12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV and Covar,"
cycle = "Light"
light_pt3.VO2.rcf <- rcf_pipe(light_pt3, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.1282987 3.2196870
##
## $cut
## [1] 0.06182777
## [1] "SPF EC50: 1.69008429724231"
## [1] "GF EC50: 1.45382056067872"
t <- light_pt3.VO2.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt3,covar_GF_pt3)) %>% set_rownames(NULL)
ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50", pt, sep = " "))
dixon.test(t[[1]])
##
## Dixon test for outliers
##
## data: t[[1]]
## Q.cohort2.3_VO2_M_9 = 0.20485, p-value = 0.9297
## alternative hypothesis: lowest value 1.43304699918151 is an outlier
dixon.test(t[[2]])
##
## Dixon test for outliers
##
## data: t[[2]]
## Q.cohort2.3_VO2_M_3 = 0.11382, p-value = 0.5661
## alternative hypothesis: highest value 1.53061348786286 is an outlier
f <- ggarrange(
ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))
shapiro.test(t[[1]])
##
## Shapiro-Wilk normality test
##
## data: t[[1]]
## W = 0.97217, p-value = 0.9126
shapiro.test(t[[2]])
##
## Shapiro-Wilk normality test
##
## data: t[[2]]
## W = 0.93438, p-value = 0.6144
var.test(t[[1]], t[[2]])
##
## F test to compare two variances
##
## data: t[[1]] and t[[2]]
## F = 6.0985, num df = 8, denom df = 5, p-value = 0.06199
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.9025158 29.3779346
## sample estimates:
## ratio of variances
## 6.098454
t.test(t[[1]], t[[2]], var.equal = T)
##
## Two Sample t-test
##
## data: t[[1]] and t[[2]]
## t = 3.5555, df = 13, p-value = 0.003519
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.08659019 0.35475341
## sample estimates:
## mean of x mean of y
## 1.677844 1.457173
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.666 -2.449 -1.202 2.489 6.212
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.012 8.379 3.343 0.00529 **
## df[[2]] 2.844 5.247 0.542 0.59695
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.129 on 13 degrees of freedom
## Multiple R-squared: 0.0221, Adjusted R-squared: -0.05312
## F-statistic: 0.2938 on 1 and 13 DF, p-value: 0.597
##
## y = 28 + 2.8x, r^2 = 0.0221
##
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 6.000000000 0.000000000 0.000000000 1.387532895 1.530613488 0.143080593
## sum median mean SE.mean CI.mean.0.95 var
## 8.743035582 1.450497011 1.457172597 0.023634555 0.060754559 0.003351553
## std.dev coef.var
## 0.057892601 0.039729405
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 1.43304700 1.89381736 0.46077036
## sum median mean SE.mean CI.mean.0.95 var
## 15.10059959 1.71904878 1.67784440 0.04765535 0.10989344 0.02043929
## std.dev coef.var
## 0.14296606 0.08520818
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 6.00000000 0.00000000 0.00000000 28.70000000 33.20000000 4.50000000
## sum median mean SE.mean CI.mean.0.95 var
## 179.60000000 29.35000000 29.93333333 0.67806915 1.74303225 2.75866667
## std.dev coef.var
## 1.66092344 0.05548742
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 30.90000000 38.30000000 7.40000000
## sum median mean SE.mean CI.mean.0.95 var
## 308.40000000 34.90000000 34.26666667 0.82276634 1.89730257 6.09250000
## std.dev coef.var
## 2.46829901 0.07203207
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.9704 0.1838
## 13
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.1753 0.17531 12.64 0.00352 **
## Residuals 13 0.1803 0.01387
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 67.60 67.60 14.05 0.00243 **
## Residuals 13 62.53 4.81
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0000011 1.100e-06 0.029 0.868
## Residuals 13 0.0004959 3.814e-05
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.17531 0.17531 24.53 0.000335 ***
## df[[3]] 1 0.09452 0.09452 13.23 0.003407 **
## Residuals 12 0.08575 0.00715
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.47390 1 66.318 3.132e-06 ***
## df[[1]] 0.26197 1 36.660 5.717e-05 ***
## df[[3]] 0.09452 1 13.227 0.003407 **
## Residuals 0.08575 12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV and Covar,"
cycle = "Dark"
dark_pt3.VO2.rcf <- rcf_pipe(dark_pt3, cycle, channel)
## [1] "Range:"
## $range
## [1] 1.133380 3.630553
##
## $cut
## [1] 0.04994346
## [1] "SPF EC50: 2.05497861966052"
## [1] "GF EC50: 1.93321699617181"
t <- dark_pt3.VO2.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt3,covar_GF_pt3)) %>% set_rownames(NULL)
ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50", pt, sep = " "))
dixon.test(t[[1]])
##
## Dixon test for outliers
##
## data: t[[1]]
## Q.cohort2.3_VO2_M_12 = 0.00219, p-value = 0.01876
## alternative hypothesis: highest value 2.35630507609501 is an outlier
dixon.test(t[[2]])
##
## Dixon test for outliers
##
## data: t[[2]]
## Q.cohort1.3_VO2_M_10 = 0.2955, p-value = 0.6757
## alternative hypothesis: lowest value 1.76718655653742 is an outlier
f <- ggarrange(
ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))
shapiro.test(t[[1]])
##
## Shapiro-Wilk normality test
##
## data: t[[1]]
## W = 0.90042, p-value = 0.2544
shapiro.test(t[[2]])
##
## Shapiro-Wilk normality test
##
## data: t[[2]]
## W = 0.91308, p-value = 0.457
var.test(t[[1]], t[[2]])
##
## F test to compare two variances
##
## data: t[[1]] and t[[2]]
## F = 3.0977, num df = 8, denom df = 5, p-value = 0.2285
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.4584277 14.9223546
## sample estimates:
## ratio of variances
## 3.097675
t.test(t[[1]], t[[2]], var.equal = T)
##
## Two Sample t-test
##
## data: t[[1]] and t[[2]]
## t = 2.1269, df = 13, p-value = 0.05315
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.002649689 0.339813505
## sample estimates:
## mean of x mean of y
## 2.088823 1.920241
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.7689 -2.7219 -0.8649 2.5665 5.8477
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.2124 10.1851 3.064 0.00904 **
## df[[2]] 0.6535 5.0225 0.130 0.89847
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.162 on 13 degrees of freedom
## Multiple R-squared: 0.001301, Adjusted R-squared: -0.07552
## F-statistic: 0.01693 on 1 and 13 DF, p-value: 0.8985
##
## y = 31 + 0.65x, r^2 = 0.0013
##
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 6.000000000 0.000000000 0.000000000 1.767186557 2.016871555 0.249684998
## sum median mean SE.mean CI.mean.0.95 var
## 11.521444188 1.941350371 1.920240698 0.040563066 0.104270682 0.009872174
## std.dev coef.var
## 0.099358815 0.051742896
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 1.87848131 2.35630508 0.47782377
## sum median mean SE.mean CI.mean.0.95 var
## 18.79940346 2.06479867 2.08882261 0.05829121 0.13441978 0.03058079
## std.dev coef.var
## 0.17487364 0.08371876
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 6.00000000 0.00000000 0.00000000 28.70000000 33.20000000 4.50000000
## sum median mean SE.mean CI.mean.0.95 var
## 179.60000000 29.35000000 29.93333333 0.67806915 1.74303225 2.75866667
## std.dev coef.var
## 1.66092344 0.05548742
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 30.90000000 38.30000000 7.40000000
## sum median mean SE.mean CI.mean.0.95 var
## 308.40000000 34.90000000 34.26666667 0.82276634 1.89730257 6.09250000
## std.dev coef.var
## 2.46829901 0.07203207
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.0775 0.3182
## 13
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.1023 0.10231 4.524 0.0531 .
## Residuals 13 0.2940 0.02262
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 67.60 67.60 14.05 0.00243 **
## Residuals 13 62.53 4.81
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0000271 2.706e-05 0.471 0.504
## Residuals 13 0.0007465 5.742e-05
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.10231 0.10231 6.015 0.0305 *
## df[[3]] 1 0.08989 0.08989 5.285 0.0403 *
## Residuals 12 0.20412 0.01701
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.64393 1 37.8564 4.926e-05 ***
## df[[1]] 0.19169 1 11.2691 0.005706 **
## df[[3]] 0.08989 1 5.2846 0.040283 *
## Residuals 0.20412 12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV and Covar,"
RQ
# RQ ####
channel = "RQ"
total_run %<>% mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. < 0.5, NA, .)) %>%
mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. > 1.5, NA, .))
total2_pt2 %<>% mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. < 0.5, NA, .)) %>%
mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. > 1.1, NA, .))
total2_pt3 %<>% mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. < 0.5, NA, .))
total_pt2 %<>% mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. < 0.5, NA, .)) %>%
mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. > 1.1, NA, .))
total_pt3 %<>% mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. < 0.5, NA, .))
light_pt2 %<>% mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. < 0.5, NA, .)) %>%
mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. > 1.1, NA, .))
dark_pt2 %<>% mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. < 0.5, NA, .)) %>%
mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. > 1.1, NA, .))
light_pt3 %<>% mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. < 0.5, NA, .))
dark_pt3 %<>% mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. < 0.5, NA, .))
p <- graph_groups(total_run, channel, title = paste0(channel, " before and after adding HF diet"))
p + geom_rect(aes(xmin=cuts[6], xmax=cuts[7], ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_rect(aes(xmin=cuts[8], xmax=cuts[9], ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_rect(aes(xmin=cuts[10], xmax=cuts[11], ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_rect(aes(xmin=cuts[12], xmax=cuts[13], ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_rect(aes(xmin=cuts[14], xmax=cuts[15], ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_line(size=1) +
geom_vline(xintercept = 1500)
# p1 <- graph_groups(total2_pt1, channel, title = paste0(channel, " before HF diet"))
# p2 <- graph_groups(total2_pt2, channel, title = paste0(channel, " immediately after HF diet"))
# p3 <- graph_groups(total2_pt3, channel, title = paste0(channel, " 1 month after HF diet"))
p1 <- graph_groups(total2_pt1, channel)
p2 <- graph_groups(total2_pt2, channel)
p3 <- graph_groups(total2_pt3, channel)
leg <- get_legend(p1)
p1 <- p1 + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
p2 <- p2 + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
p3 <- p3 + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
f <- ggarrange(p1, p2, p3, leg)
annotate_figure(f, top = text_grob(paste0(channel, ", before, after, and a month after adding HF diet"), face = "bold"), left = text_grob(channel), bottom = text_grob("observation"))
Before HF diet
# pt1 RQ ####
covar_SPF_pt1 <- c(cohort1_1_covar[[1]][c(2,3)], cohort2_1_covar[[1]][c(6,7,8,9)], cohort3_1_covar[[1]][c(4,5,6)])
covar_GF_pt1 <- c(cohort1_1_covar[[1]][1], cohort2_1_covar[[1]][c(1,2,3,4,5)], cohort3_1_covar[[1]][c(1,2,3)])
pt = "before HF diet"
cycle = "Total"
total_pt1.RQ.rcf <- rcf_pipe(total_pt1, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.2911567 1.3067190
##
## $cut
## [1] 0.02031125
## [1] "SPF EC50: 0.86888794275416"
## [1] "GF EC50: 0.826388799709833"
t <- total_pt1.RQ.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)
ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50", pt, sep = " "))
dixon.test(t[[1]])
##
## Dixon test for outliers
##
## data: t[[1]]
## Q.cohort2.1_RQ_M_9 = 0.48246, p-value = 0.1359
## alternative hypothesis: lowest value 0.823158188054464 is an outlier
dixon.test(t[[2]])
##
## Dixon test for outliers
##
## data: t[[2]]
## Q.cohort3.1_RQ_M_9 = 0.56455, p-value = 0.0536
## alternative hypothesis: lowest value 0.760142288370411 is an outlier
f <- ggarrange(
ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", sep = " "), face = "bold"))
shapiro.test(t[[1]])
##
## Shapiro-Wilk normality test
##
## data: t[[1]]
## W = 0.95552, p-value = 0.7505
shapiro.test(t[[2]])
##
## Shapiro-Wilk normality test
##
## data: t[[2]]
## W = 0.86605, p-value = 0.1115
var.test(t[[1]], t[[2]])
##
## F test to compare two variances
##
## data: t[[1]] and t[[2]]
## F = 0.73521, num df = 8, denom df = 8, p-value = 0.6738
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.1658385 3.2593552
## sample estimates:
## ratio of variances
## 0.7352051
# t.test(t[[1]], t[[2]])
# If var.test is not significant, set var.equal to true!
t.test(t[[1]], t[[2]], var.equal = T)
##
## Two Sample t-test
##
## data: t[[1]] and t[[2]]
## t = 3.2744, df = 16, p-value = 0.00477
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.01572764 0.07348722
## sample estimates:
## mean of x mean of y
## 0.8719954 0.8273880
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4906 -1.1167 -0.3626 1.0107 5.2180
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 51.85 11.83 4.382 0.000465 ***
## df[[2]] -26.46 13.91 -1.902 0.075384 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.079 on 16 degrees of freedom
## Multiple R-squared: 0.1843, Adjusted R-squared: 0.1334
## F-statistic: 3.616 on 1 and 16 DF, p-value: 0.07538
##
## y = 52 + -26x, r^2 = 0.184
##
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 9.0000000000 0.0000000000 0.0000000000 0.7601422884 0.8762337005 0.1160914122
## sum median mean SE.mean CI.mean.0.95 var
## 7.4464918248 0.8317564691 0.8273879805 0.0103419439 0.0238485654 0.0009626022
## std.dev coef.var
## 0.0310258317 0.0374985284
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.000000000 0.000000000 0.000000000 0.823158188 0.914947835 0.091789646
## sum median mean SE.mean CI.mean.0.95 var
## 7.847958723 0.864325085 0.871995414 0.008867607 0.020448738 0.000707710
## std.dev coef.var
## 0.026602820 0.030507982
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 29.00000000 35.00000000 6.00000000
## sum median mean SE.mean CI.mean.0.95 var
## 276.00000000 30.00000000 30.66666667 0.68211273 1.57295478 4.18750000
## std.dev coef.var
## 2.04633819 0.06672842
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 25.60000000 30.40000000 4.80000000
## sum median mean SE.mean CI.mean.0.95 var
## 252.60000000 27.80000000 28.06666667 0.53800041 1.24063118 2.60500000
## std.dev coef.var
## 1.61400124 0.05750598
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.0132 0.91
## 16
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.008954 0.008954 10.72 0.00477 **
## Residuals 16 0.013362 0.000835
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 30.42 30.420 8.957 0.00861 **
## Residuals 16 54.34 3.396
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 7.502e-05 7.502e-05 16.98 0.000801 ***
## Residuals 16 7.069e-05 4.420e-06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.008954 0.008954 10.117 0.0062 **
## df[[3]] 1 0.000087 0.000087 0.098 0.7587
## Residuals 15 0.013276 0.000885
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.043067 1 48.6603 4.458e-06 ***
## df[[1]] 0.004927 1 5.5667 0.03228 *
## df[[3]] 0.000087 1 0.0979 0.75872
## Residuals 0.013276 15
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.0322840364003511"
cycle = "Light"
light_pt1.RQ.rcf <- rcf_pipe(light_pt1, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.4186665 0.9756889
##
## $cut
## [1] 0.01114045
## [1] "SPF EC50: 0.808304086759372"
## [1] "GF EC50: 0.756862038448115"
t <- light_pt1.RQ.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)
ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50", pt, sep = " "))
dixon.test(t[[1]])
##
## Dixon test for outliers
##
## data: t[[1]]
## Q.cohort2.1_RQ_M_11 = 0.12043, p-value = 0.6669
## alternative hypothesis: highest value 0.855640198857886 is an outlier
dixon.test(t[[2]])
##
## Dixon test for outliers
##
## data: t[[2]]
## Q.cohort3.1_RQ_M_9 = 0.45221, p-value = 0.181
## alternative hypothesis: lowest value 0.70776270316176 is an outlier
f <- ggarrange(
ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))
shapiro.test(t[[1]])
##
## Shapiro-Wilk normality test
##
## data: t[[1]]
## W = 0.89747, p-value = 0.2377
shapiro.test(t[[2]])
##
## Shapiro-Wilk normality test
##
## data: t[[2]]
## W = 0.91972, p-value = 0.3899
var.test(t[[1]], t[[2]])
##
## F test to compare two variances
##
## data: t[[1]] and t[[2]]
## F = 2.4598, num df = 8, denom df = 8, p-value = 0.2245
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.5548409 10.9047281
## sample estimates:
## ratio of variances
## 2.459754
t.test(t[[1]], t[[2]], var.equal = T)
##
## Two Sample t-test
##
## data: t[[1]] and t[[2]]
## t = 3.2633, df = 16, p-value = 0.004882
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.01710230 0.08051796
## sample estimates:
## mean of x mean of y
## 0.8054247 0.7566145
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.1339 -0.8929 -0.3644 0.7451 4.9434
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 47.952 9.955 4.817 0.00019 ***
## df[[2]] -23.796 12.731 -1.869 0.08001 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.085 on 16 degrees of freedom
## Multiple R-squared: 0.1792, Adjusted R-squared: 0.1279
## F-statistic: 3.494 on 1 and 16 DF, p-value: 0.08001
##
## y = 48 + -24x, r^2 = 0.179
##
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 9.0000000000 0.0000000000 0.0000000000 0.7077627032 0.7865111501 0.0787484469
## sum median mean SE.mean CI.mean.0.95 var
## 6.8095307367 0.7534788906 0.7566145263 0.0080413220 0.0185433218 0.0005819657
## std.dev coef.var
## 0.0241239660 0.0318840904
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.000000000 0.000000000 0.000000000 0.756513163 0.855640199 0.099127036
## sum median mean SE.mean CI.mean.0.95 var
## 7.248821878 0.814399686 0.805424653 0.012611690 0.029082609 0.001431492
## std.dev coef.var
## 0.037835069 0.046975306
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 29.00000000 35.00000000 6.00000000
## sum median mean SE.mean CI.mean.0.95 var
## 276.00000000 30.00000000 30.66666667 0.68211273 1.57295478 4.18750000
## std.dev coef.var
## 2.04633819 0.06672842
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 25.60000000 30.40000000 4.80000000
## sum median mean SE.mean CI.mean.0.95 var
## 252.60000000 27.80000000 28.06666667 0.53800041 1.24063118 2.60500000
## std.dev coef.var
## 1.61400124 0.05750598
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 3.5522 0.07776 .
## 16
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.01072 0.010721 10.65 0.00488 **
## Residuals 16 0.01611 0.001007
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 30.42 30.420 8.957 0.00861 **
## Residuals 16 54.34 3.396
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 7.257e-05 7.257e-05 17.99 0.000622 ***
## Residuals 16 6.454e-05 4.030e-06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.010721 0.010721 10.036 0.00637 **
## df[[3]] 1 0.000083 0.000083 0.078 0.78371
## Residuals 15 0.016024 0.001068
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.036251 1 33.9339 3.342e-05 ***
## df[[1]] 0.005996 1 5.6126 0.03168 *
## df[[3]] 0.000083 1 0.0781 0.78371
## Residuals 0.016024 15
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.0316788344413115"
cycle = "Dark"
dark_pt1.RQ.rcf <- rcf_pipe(dark_pt1, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.2911567 1.3067190
##
## $cut
## [1] 0.02031125
## [1] "SPF EC50: 0.905887016006714"
## [1] "GF EC50: 0.888108777077216"
t <- dark_pt1.RQ.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)
ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50", pt, sep = " "))
dixon.test(t[[1]])
##
## Dixon test for outliers
##
## data: t[[1]]
## Q.cohort3.1_RQ_M_12 = 0.44552, p-value = 0.1922
## alternative hypothesis: lowest value 0.853686998799925 is an outlier
dixon.test(t[[2]])
##
## Dixon test for outliers
##
## data: t[[2]]
## Q.cohort3.1_RQ_M_9 = 0.19108, p-value = 0.9903
## alternative hypothesis: lowest value 0.820799851235852 is an outlier
f <- ggarrange(
ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))
shapiro.test(t[[1]])
##
## Shapiro-Wilk normality test
##
## data: t[[1]]
## W = 0.96529, p-value = 0.8518
shapiro.test(t[[2]])
##
## Shapiro-Wilk normality test
##
## data: t[[2]]
## W = 0.90587, p-value = 0.288
var.test(t[[1]], t[[2]])
##
## F test to compare two variances
##
## data: t[[1]] and t[[2]]
## F = 0.68623, num df = 8, denom df = 8, p-value = 0.6068
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.1547903 3.0422174
## sample estimates:
## ratio of variances
## 0.6862258
t.test(t[[1]], t[[2]], var.equal = T)
##
## Two Sample t-test
##
## data: t[[1]] and t[[2]]
## t = 1.8123, df = 16, p-value = 0.08875
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.004346843 0.055560099
## sample estimates:
## mean of x mean of y
## 0.9066351 0.8810285
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6795 -1.2565 -0.2520 0.8756 5.8246
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.77 15.33 2.725 0.015 *
## df[[2]] -13.87 17.14 -0.809 0.430
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.256 on 16 degrees of freedom
## Multiple R-squared: 0.03933, Adjusted R-squared: -0.02071
## F-statistic: 0.6551 on 1 and 16 DF, p-value: 0.4302
##
## y = 42 + -14x, r^2 = 0.0393
##
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 9.000000000 0.000000000 0.000000000 0.820799851 0.919892794 0.099092943
## sum median mean SE.mean CI.mean.0.95 var
## 7.929256641 0.890881776 0.881028516 0.010881101 0.025091863 0.001065585
## std.dev coef.var
## 0.032643302 0.037051357
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.0000000000 0.0000000000 0.0000000000 0.8536869988 0.9504519265 0.0967649277
## sum median mean SE.mean CI.mean.0.95 var
## 8.1597162900 0.9084909024 0.9066351433 0.0090137676 0.0207857854 0.0007312321
## std.dev coef.var
## 0.0270413029 0.0298260034
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 29.00000000 35.00000000 6.00000000
## sum median mean SE.mean CI.mean.0.95 var
## 276.00000000 30.00000000 30.66666667 0.68211273 1.57295478 4.18750000
## std.dev coef.var
## 2.04633819 0.06672842
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 25.60000000 30.40000000 4.80000000
## sum median mean SE.mean CI.mean.0.95 var
## 252.60000000 27.80000000 28.06666667 0.53800041 1.24063118 2.60500000
## std.dev coef.var
## 1.61400124 0.05750598
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.196 0.6639
## 16
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.002951 0.0029506 3.284 0.0888 .
## Residuals 16 0.014375 0.0008984
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 30.42 30.420 8.957 0.00861 **
## Residuals 16 54.34 3.396
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 5.744e-05 5.744e-05 13.1 0.0023 **
## Residuals 16 7.016e-05 4.390e-06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.002951 0.0029506 3.093 0.099 .
## df[[3]] 1 0.000065 0.0000646 0.068 0.798
## Residuals 15 0.014310 0.0009540
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.041245 1 43.2339 8.805e-06 ***
## df[[1]] 0.002334 1 2.4464 0.1386
## df[[3]] 0.000065 1 0.0678 0.7982
## Residuals 0.014310 15
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
Immediately after HF diet
# pt2 RQ ####
# delete these once I add cohort1_2!!
covar_SPF_pt1 <- c(cohort2_1_covar[[1]][c(6,7,8,9)], cohort3_1_covar[[1]][c(4,5,6)])
covar_GF_pt1 <- c(cohort2_1_covar[[1]][c(1,2,3,4,5)], cohort3_1_covar[[1]][c(1,2,3)])
pt = "immediately after HF diet"
cycle = "Total"
total_pt2.RQ.rcf <- rcf_pipe(total_pt2, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.5205333 1.0818660
##
## $cut
## [1] 0.01122665
## [1] "SPF EC50: 0.889997042895635"
## [1] "GF EC50: 0.813386849541187"
t <- total_pt2.RQ.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)
ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50", pt, sep = " "))
dixon.test(t[[1]])
##
## Dixon test for outliers
##
## data: t[[1]]
## Q.cohort2.2_RQ_M_9 = 0.4003, p-value = 0.2641
## alternative hypothesis: lowest value 0.841189809990053 is an outlier
dixon.test(t[[2]])
##
## Dixon test for outliers
##
## data: t[[2]]
## Q.cohort3.2_RQ_M_9 = 0.20498, p-value = 0.9789
## alternative hypothesis: lowest value 0.724630874890402 is an outlier
f <- ggarrange(
ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))
shapiro.test(t[[1]])
##
## Shapiro-Wilk normality test
##
## data: t[[1]]
## W = 0.90654, p-value = 0.3725
shapiro.test(t[[2]])
##
## Shapiro-Wilk normality test
##
## data: t[[2]]
## W = 0.91694, p-value = 0.4055
var.test(t[[1]], t[[2]])
##
## F test to compare two variances
##
## data: t[[1]] and t[[2]]
## F = 0.25911, num df = 6, denom df = 7, p-value = 0.1205
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.0506214 1.4757565
## sample estimates:
## ratio of variances
## 0.2591105
t.test(t[[1]], t[[2]], var.equal = FALSE) # F-test p-val: 0.03675
##
## Welch Two Sample t-test
##
## data: t[[1]] and t[[2]]
## t = 4.0728, df = 10.668, p-value = 0.001961
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.0358323 0.1208002
## sample estimates:
## mean of x mean of y
## 0.8851837 0.8068675
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5710 -1.2369 -0.4368 1.0528 4.4375
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 44.242 8.279 5.344 0.000133 ***
## df[[2]] -17.210 9.796 -1.757 0.102472
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.018 on 13 degrees of freedom
## Multiple R-squared: 0.1919, Adjusted R-squared: 0.1297
## F-statistic: 3.086 on 1 and 13 DF, p-value: 0.1025
##
## y = 44 + -17x, r^2 = 0.192
##
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 8.000000000 0.000000000 0.000000000 0.724630875 0.860256238 0.135625364
## sum median mean SE.mean CI.mean.0.95 var
## 6.454939799 0.819557590 0.806867475 0.016890152 0.039938864 0.002282218
## std.dev coef.var
## 0.047772565 0.059207449
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 7.0000000000 0.0000000000 0.0000000000 0.8411898100 0.9090888065 0.0678989965
## sum median mean SE.mean CI.mean.0.95 var
## 6.1962861453 0.8931365901 0.8851837350 0.0091911972 0.0224900493 0.0005913467
## std.dev coef.var
## 0.0243176220 0.0274718355
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 8.00000000 0.00000000 0.00000000 29.00000000 35.00000000 6.00000000
## sum median mean SE.mean CI.mean.0.95 var
## 246.50000000 30.00000000 30.81250000 0.75555975 1.78661491 4.56696429
## std.dev coef.var
## 2.13704569 0.06935645
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 7.0000000 0.0000000 0.0000000 26.3000000 30.4000000 4.1000000
## sum median mean SE.mean CI.mean.0.95 var
## 199.4000000 28.1000000 28.4857143 0.5629127 1.3773978 2.2180952
## std.dev coef.var
## 1.4893271 0.0522833
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.4961 0.243
## 13
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.02290 0.022898 15.25 0.00181 **
## Residuals 13 0.01952 0.001502
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 20.21 20.212 5.803 0.0315 *
## Residuals 13 45.28 3.483
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 8.775e-05 8.775e-05 18.31 0.000897 ***
## Residuals 13 6.230e-05 4.790e-06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.022898 0.022898 14.114 0.00274 **
## df[[3]] 1 0.000055 0.000055 0.034 0.85740
## Residuals 12 0.019469 0.001622
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.033509 1 20.6538 0.0006724 ***
## df[[1]] 0.014814 1 9.1309 0.0106284 *
## df[[3]] 0.000055 1 0.0337 0.8573954
## Residuals 0.019469 12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.0106283891377386"
cycle = "Light"
light_pt2.RQ.rcf <- rcf_pipe(light_pt2, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.5205333 1.0163130
##
## $cut
## [1] 0.009915594
## [1] "SPF EC50: 0.866984758589404"
## [1] "GF EC50: 0.787869144945168"
t <- light_pt2.RQ.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)
ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50", pt, sep = " "))
dixon.test(t[[1]])
##
## Dixon test for outliers
##
## data: t[[1]]
## Q.cohort3.2_RQ_M_12 = 0.31499, p-value = 0.4858
## alternative hypothesis: lowest value 0.810903257955022 is an outlier
dixon.test(t[[2]])
##
## Dixon test for outliers
##
## data: t[[2]]
## Q.cohort3.2_RQ_M_9 = 0.39154, p-value = 0.3827
## alternative hypothesis: lowest value 0.70694525860472 is an outlier
f <- ggarrange(
ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))
shapiro.test(t[[1]])
##
## Shapiro-Wilk normality test
##
## data: t[[1]]
## W = 0.95353, p-value = 0.7617
shapiro.test(t[[2]])
##
## Shapiro-Wilk normality test
##
## data: t[[2]]
## W = 0.9858, p-value = 0.9858
var.test(t[[1]], t[[2]])
##
## F test to compare two variances
##
## data: t[[1]] and t[[2]]
## F = 0.5345, num df = 6, denom df = 7, p-value = 0.463
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.1044234 3.0442349
## sample estimates:
## ratio of variances
## 0.534501
t.test(t[[1]], t[[2]], var.equal = T)
##
## Two Sample t-test
##
## data: t[[1]] and t[[2]]
## t = 3.5647, df = 13, p-value = 0.003457
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.03017021 0.12299306
## sample estimates:
## mean of x mean of y
## 0.8638716 0.7872899
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5778 -1.2277 -0.5805 1.1171 4.2483
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 45.344 7.651 5.927 5.01e-05 ***
## df[[2]] -18.976 9.275 -2.046 0.0616 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.952 on 13 degrees of freedom
## Multiple R-squared: 0.2435, Adjusted R-squared: 0.1854
## F-statistic: 4.185 on 1 and 13 DF, p-value: 0.06156
##
## y = 45 + -19x, r^2 = 0.244
##
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 8.000000000 0.000000000 0.000000000 0.706945259 0.860038846 0.153093587
## sum median mean SE.mean CI.mean.0.95 var
## 6.298319331 0.787148223 0.787289916 0.016562386 0.039163820 0.002194501
## std.dev coef.var
## 0.046845502 0.059502225
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 7.000000000 0.000000000 0.000000000 0.810903258 0.905564075 0.094660817
## sum median mean SE.mean CI.mean.0.95 var
## 6.047100861 0.867761027 0.863871552 0.012944735 0.031674625 0.001172963
## std.dev coef.var
## 0.034248549 0.039645418
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 8.00000000 0.00000000 0.00000000 29.00000000 35.00000000 6.00000000
## sum median mean SE.mean CI.mean.0.95 var
## 246.50000000 30.00000000 30.81250000 0.75555975 1.78661491 4.56696429
## std.dev coef.var
## 2.13704569 0.06935645
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 7.0000000 0.0000000 0.0000000 26.3000000 30.4000000 4.1000000
## sum median mean SE.mean CI.mean.0.95 var
## 199.4000000 28.1000000 28.4857143 0.5629127 1.3773978 2.2180952
## std.dev coef.var
## 1.4893271 0.0522833
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.3328 0.5739
## 13
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0219 0.021895 12.71 0.00346 **
## Residuals 13 0.0224 0.001723
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 20.21 20.212 5.803 0.0315 *
## Residuals 13 45.28 3.483
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 8.345e-05 8.345e-05 15.39 0.00175 **
## Residuals 13 7.051e-05 5.420e-06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.021895 0.021895 12.096 0.00456 **
## df[[3]] 1 0.000679 0.000679 0.375 0.55179
## Residuals 12 0.021721 0.001810
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.038963 1 21.5255 0.0005706 ***
## df[[1]] 0.011786 1 6.5114 0.0253839 *
## df[[3]] 0.000679 1 0.3749 0.5517931
## Residuals 0.021721 12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.0253839025719477"
# cohort2.2_RQ_M_9 is probably an outlier
# dark_pt2 %<>% select(-cohort2.2_RQ_M_9)
# covar_SPF_pt1 <- c(cohort2_1_covar[[1]][c(6,8,9)], cohort3_1_covar[[1]][c(4,5,6)])
cycle = "Dark"
dark_pt2.RQ.rcf <- rcf_pipe(dark_pt2, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.6168776 1.0818660
##
## $cut
## [1] 0.009299768
## [1] "SPF EC50: 0.902146907095478"
## [1] "GF EC50: 0.832689248018978"
t <- dark_pt2.RQ.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)
ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50", pt, sep = " "))
dixon.test(t[[1]])
##
## Dixon test for outliers
##
## data: t[[1]]
## Q.cohort2.2_RQ_M_9 = 0.67752, p-value = 0.01046
## alternative hypothesis: lowest value 0.84275526491303 is an outlier
dixon.test(t[[2]])
##
## Dixon test for outliers
##
## data: t[[2]]
## Q.cohort3.2_RQ_M_9 = 0.073739, p-value = 0.3761
## alternative hypothesis: lowest value 0.744699285944684 is an outlier
f <- ggarrange(
ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))
shapiro.test(t[[1]])
##
## Shapiro-Wilk normality test
##
## data: t[[1]]
## W = 0.74009, p-value = 0.009986
shapiro.test(t[[2]])
##
## Shapiro-Wilk normality test
##
## data: t[[2]]
## W = 0.85735, p-value = 0.113
var.test(t[[1]], t[[2]])
##
## F test to compare two variances
##
## data: t[[1]] and t[[2]]
## F = 0.2754, num df = 6, denom df = 7, p-value = 0.1371
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.05380359 1.56852627
## sample estimates:
## ratio of variances
## 0.2753989
t.test(t[[1]], t[[2]], var.equal = T)
##
## Two Sample t-test
##
## data: t[[1]] and t[[2]]
## t = 3.6963, df = 13, p-value = 0.002689
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.0313418 0.1195076
## sample estimates:
## mean of x mean of y
## 0.8969731 0.8215484
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7522 -1.2915 -0.6480 0.9957 4.7798
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 42.167 8.809 4.787 0.000355 ***
## df[[2]] -14.521 10.263 -1.415 0.180605
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.089 on 13 degrees of freedom
## Multiple R-squared: 0.1334, Adjusted R-squared: 0.06679
## F-statistic: 2.002 on 1 and 13 DF, p-value: 0.1806
##
## y = 42 + -15x, r^2 = 0.133
##
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 8.000000000 0.000000000 0.000000000 0.744699286 0.871081793 0.126382507
## sum median mean SE.mean CI.mean.0.95 var
## 6.572387075 0.834724603 0.821548384 0.017086348 0.040402793 0.002335546
## std.dev coef.var
## 0.048327490 0.058824886
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 7.0000000000 0.0000000000 0.0000000000 0.8427552649 0.9241419801 0.0813867152
## sum median mean SE.mean CI.mean.0.95 var
## 6.2788116371 0.9031934811 0.8969730910 0.0095857548 0.0234554971 0.0006432069
## std.dev coef.var
## 0.0253615234 0.0282745644
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 8.00000000 0.00000000 0.00000000 29.00000000 35.00000000 6.00000000
## sum median mean SE.mean CI.mean.0.95 var
## 246.50000000 30.00000000 30.81250000 0.75555975 1.78661491 4.56696429
## std.dev coef.var
## 2.13704569 0.06935645
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 7.0000000 0.0000000 0.0000000 26.3000000 30.4000000 4.1000000
## sum median mean SE.mean CI.mean.0.95 var
## 199.4000000 28.1000000 28.4857143 0.5629127 1.3773978 2.2180952
## std.dev coef.var
## 1.4893271 0.0522833
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 2.815 0.1172
## 13
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.02124 0.021239 13.66 0.00269 **
## Residuals 13 0.02021 0.001554
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 20.21 20.212 5.803 0.0315 *
## Residuals 13 45.28 3.483
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 8.558e-05 8.558e-05 18.97 0.000779 ***
## Residuals 13 5.865e-05 4.510e-06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.021239 0.021239 12.651 0.00395 **
## df[[3]] 1 0.000063 0.000063 0.037 0.84979
## Residuals 12 0.020145 0.001679
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.0292315 1 17.4125 0.001293 **
## df[[1]] 0.0157706 1 9.3941 0.009808 **
## df[[3]] 0.0000629 1 0.0374 0.849794
## Residuals 0.0201452 12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.00980775246066405"
1 month after HF diet
# pt3 RQ ####
# cohort2.3_RQ_M_1 is an outlier.
total_pt3 %<>% select(-cohort2.3_RQ_M_1)
light_pt3 %<>% select(-cohort2.3_RQ_M_1)
dark_pt3 %<>% select(-cohort2.3_RQ_M_1)
covar_SPF_pt3 <- c(cohort1_3_covar[[1]][c(2,3)], cohort2_3_covar[[1]][c(4,5,6,7)], cohort3_3_covar[[1]][c(4,5,6)])
covar_GF_pt3 <- c(cohort1_3_covar[[1]][1], cohort2_3_covar[[1]][c(2,3)], cohort3_3_covar[[1]][c(1,2,3)])
pt = "1 month after HF diet"
cycle = "Total"
total_pt3.RQ.rcf <- rcf_pipe(total_pt3, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.5034941 1.3791840
##
## $cut
## [1] 0.0175138
## [1] "SPF EC50: 0.843815968423415"
## [1] "GF EC50: 0.814793610952982"
t <- total_pt3.RQ.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt3,covar_GF_pt3)) %>% set_rownames(NULL)
ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50", pt, sep = " "))
dixon.test(t[[1]])
##
## Dixon test for outliers
##
## data: t[[1]]
## Q.cohort2.3_RQ_M_12 = 0.3487, p-value = 0.4091
## alternative hypothesis: lowest value 0.785868886334165 is an outlier
dixon.test(t[[2]])
##
## Dixon test for outliers
##
## data: t[[2]]
## Q.cohort3.3_RQ_M_2 = 0.55225, p-value = 0.1078
## alternative hypothesis: lowest value 0.779598208838337 is an outlier
f <- ggarrange(
ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))
shapiro.test(t[[1]])
##
## Shapiro-Wilk normality test
##
## data: t[[1]]
## W = 0.9508, p-value = 0.6989
shapiro.test(t[[2]])
##
## Shapiro-Wilk normality test
##
## data: t[[2]]
## W = 0.90215, p-value = 0.3868
var.test(t[[1]], t[[2]])
##
## F test to compare two variances
##
## data: t[[1]] and t[[2]]
## F = 2.1499, num df = 8, denom df = 5, p-value = 0.4148
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.3181729 10.3568953
## sample estimates:
## ratio of variances
## 2.149949
t.test(t[[1]], t[[2]], var.equal = T)
##
## Two Sample t-test
##
## data: t[[1]] and t[[2]]
## t = 1.8393, df = 13, p-value = 0.08881
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.004607271 0.057401018
## sample estimates:
## mean of x mean of y
## 0.8425080 0.8161111
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.571 -2.741 -1.206 2.654 5.123
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.47 23.36 0.620 0.546
## df[[2]] 21.71 28.07 0.773 0.453
##
## Residual standard error: 3.094 on 13 degrees of freedom
## Multiple R-squared: 0.04398, Adjusted R-squared: -0.02956
## F-statistic: 0.5981 on 1 and 13 DF, p-value: 0.4531
##
## y = 14 + 22x, r^2 = 0.044
##
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 6.0000000000 0.0000000000 0.0000000000 0.7795982088 0.8372413178 0.0576431090
## sum median mean SE.mean CI.mean.0.95 var
## 4.8966667261 0.8173628859 0.8161111210 0.0085068096 0.0218674502 0.0004341949
## std.dev coef.var
## 0.0208373428 0.0255324824
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.0000000000 0.0000000000 0.0000000000 0.7858688863 0.8876941418 0.1018252554
## sum median mean SE.mean CI.mean.0.95 var
## 7.5825719507 0.8499652336 0.8425079945 0.0101843925 0.0234852513 0.0009334967
## std.dev coef.var
## 0.0305531776 0.0362645551
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 6.00000000 0.00000000 0.00000000 28.70000000 33.20000000 4.50000000
## sum median mean SE.mean CI.mean.0.95 var
## 179.60000000 29.35000000 29.93333333 0.67806915 1.74303225 2.75866667
## std.dev coef.var
## 1.66092344 0.05548742
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 30.90000000 38.30000000 7.40000000
## sum median mean SE.mean CI.mean.0.95 var
## 308.40000000 34.90000000 34.26666667 0.82276634 1.89730257 6.09250000
## std.dev coef.var
## 2.46829901 0.07203207
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.5224 0.4826
## 13
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.002508 0.0025085 3.383 0.0888 .
## Residuals 13 0.009639 0.0007415
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 67.60 67.60 14.05 0.00243 **
## Residuals 13 62.53 4.81
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 2.523e-05 2.523e-05 6.242 0.0267 *
## Residuals 13 5.254e-05 4.042e-06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.002508 0.0025085 3.241 0.097 .
## df[[3]] 1 0.000351 0.0003508 0.453 0.514
## Residuals 12 0.009288 0.0007740
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.054279 1 70.1275 2.347e-06 ***
## df[[1]] 0.002325 1 3.0038 0.1087
## df[[3]] 0.000351 1 0.4532 0.5136
## Residuals 0.009288 12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
cycle = "Light"
light_pt3.RQ.rcf <- rcf_pipe(light_pt3, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.5164798 1.3690690
##
## $cut
## [1] 0.01705178
## [1] "SPF EC50: 0.796725491719892"
## [1] "GF EC50: 0.750865595089389"
t <- light_pt3.RQ.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt3,covar_GF_pt3)) %>% set_rownames(NULL)
ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50", pt, sep = " "))
dixon.test(t[[1]])
##
## Dixon test for outliers
##
## data: t[[1]]
## Q.cohort2.3_RQ_M_12 = 0.21909, p-value = 0.8691
## alternative hypothesis: lowest value 0.733507805063572 is an outlier
dixon.test(t[[2]])
##
## Dixon test for outliers
##
## data: t[[2]]
## Q.cohort2.3_RQ_M_3 = 0.11486, p-value = 0.5711
## alternative hypothesis: highest value 0.821379469025812 is an outlier
f <- ggarrange(
ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))
shapiro.test(t[[1]])
##
## Shapiro-Wilk normality test
##
## data: t[[1]]
## W = 0.98937, p-value = 0.9952
shapiro.test(t[[2]])
##
## Shapiro-Wilk normality test
##
## data: t[[2]]
## W = 0.90861, p-value = 0.4273
var.test(t[[1]], t[[2]])
##
## F test to compare two variances
##
## data: t[[1]] and t[[2]]
## F = 0.76136, num df = 8, denom df = 5, p-value = 0.6956
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.1126738 3.6676638
## sample estimates:
## ratio of variances
## 0.7613565
t.test(t[[1]], t[[2]], var.equal = T)
##
## Two Sample t-test
##
## data: t[[1]] and t[[2]]
## t = 1.666, df = 13, p-value = 0.1196
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.01087743 0.08418140
## sample estimates:
## mean of x mean of y
## 0.7969794 0.7603274
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.423 -2.636 -1.206 2.447 5.490
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.72 14.59 1.421 0.179
## df[[2]] 15.10 18.62 0.811 0.432
##
## Residual standard error: 3.087 on 13 degrees of freedom
## Multiple R-squared: 0.04815, Adjusted R-squared: -0.02507
## F-statistic: 0.6576 on 1 and 13 DF, p-value: 0.432
##
## y = 21 + 15x, r^2 = 0.0481
##
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 6.000000000 0.000000000 0.000000000 0.702764067 0.821379469 0.118615402
## sum median mean SE.mean CI.mean.0.95 var
## 4.561964193 0.746089137 0.760327366 0.018450130 0.047427568 0.002042444
## std.dev coef.var
## 0.045193403 0.059439401
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.000000000 0.000000000 0.000000000 0.733507805 0.860417500 0.126909695
## sum median mean SE.mean CI.mean.0.95 var
## 7.172814154 0.800611733 0.796979350 0.013144613 0.030311532 0.001555028
## std.dev coef.var
## 0.039433839 0.049479123
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 6.00000000 0.00000000 0.00000000 28.70000000 33.20000000 4.50000000
## sum median mean SE.mean CI.mean.0.95 var
## 179.60000000 29.35000000 29.93333333 0.67806915 1.74303225 2.75866667
## std.dev coef.var
## 1.66092344 0.05548742
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 30.90000000 38.30000000 7.40000000
## sum median mean SE.mean CI.mean.0.95 var
## 308.40000000 34.90000000 34.26666667 0.82276634 1.89730257 6.09250000
## std.dev coef.var
## 2.46829901 0.07203207
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.0367 0.8511
## 13
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.004836 0.004836 2.775 0.12
## Residuals 13 0.022652 0.001742
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 67.60 67.60 14.05 0.00243 **
## Residuals 13 62.53 4.81
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 1.648e-05 1.648e-05 3.435 0.0866 .
## Residuals 13 6.236e-05 4.797e-06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.004836 0.004836 2.607 0.132
## df[[3]] 1 0.000393 0.000393 0.212 0.654
## Residuals 12 0.022259 0.001855
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.048143 1 25.9538 0.0002642 ***
## df[[1]] 0.003906 1 2.1055 0.1724047
## df[[3]] 0.000393 1 0.2119 0.6535429
## Residuals 0.022259 12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
cycle = "Dark"
dark_pt3.RQ.rcf <- rcf_pipe(dark_pt3, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.5034941 1.3791840
##
## $cut
## [1] 0.0175138
## [1] "SPF EC50: 0.871701952811977"
## [1] "GF EC50: 0.866975323021419"
t <- dark_pt3.RQ.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt3,covar_GF_pt3)) %>% set_rownames(NULL)
ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50", pt, sep = " "))
dixon.test(t[[1]])
##
## Dixon test for outliers
##
## data: t[[1]]
## Q.cohort2.3_RQ_M_10 = 0.041949, p-value = 0.25
## alternative hypothesis: lowest value 0.831931716444821 is an outlier
dixon.test(t[[2]])
##
## Dixon test for outliers
##
## data: t[[2]]
## Q.cohort3.3_RQ_M_1 = 0.44814, p-value = 0.2602
## alternative hypothesis: lowest value 0.805485833437532 is an outlier
f <- ggarrange(
ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))
shapiro.test(t[[1]])
##
## Shapiro-Wilk normality test
##
## data: t[[1]]
## W = 0.90685, p-value = 0.2944
shapiro.test(t[[2]])
##
## Shapiro-Wilk normality test
##
## data: t[[2]]
## W = 0.92123, p-value = 0.5142
var.test(t[[1]], t[[2]])
##
## F test to compare two variances
##
## data: t[[1]] and t[[2]]
## F = 0.53355, num df = 8, denom df = 5, p-value = 0.4096
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.07896096 2.57027073
## sample estimates:
## ratio of variances
## 0.5335528
t.test(t[[1]], t[[2]], var.equal = T)
##
## Two Sample t-test
##
## data: t[[1]] and t[[2]]
## t = 0.40859, df = 13, p-value = 0.6895
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.02685763 0.03938625
## sample estimates:
## mean of x mean of y
## 0.8704050 0.8641407
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9276 -2.7650 -0.7837 2.5419 5.8042
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.84 26.02 1.30 0.216
## df[[2]] -1.51 29.97 -0.05 0.961
##
## Residual standard error: 3.164 on 13 degrees of freedom
## Multiple R-squared: 0.0001953, Adjusted R-squared: -0.07671
## F-statistic: 0.002539 on 1 and 13 DF, p-value: 0.9606
##
## y = 34 + -1.5x, r^2 = 0.000195
##
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 6.000000000 0.000000000 0.000000000 0.805485833 0.899655802 0.094169968
## sum median mean SE.mean CI.mean.0.95 var
## 5.184844213 0.871920663 0.864140702 0.014064760 0.036154616 0.001186905
## std.dev coef.var
## 0.034451485 0.039867911
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.0000000000 0.0000000000 0.0000000000 0.8319317164 0.9007889218 0.0688572053
## sum median mean SE.mean CI.mean.0.95 var
## 7.8336450909 0.8710080145 0.8704050101 0.0083883275 0.0193435180 0.0006332764
## std.dev coef.var
## 0.0251649826 0.0289118081
## Covar - IV:
## df[[1]]: GF
## nbr.val nbr.null nbr.na min max range
## 6.00000000 0.00000000 0.00000000 28.70000000 33.20000000 4.50000000
## sum median mean SE.mean CI.mean.0.95 var
## 179.60000000 29.35000000 29.93333333 0.67806915 1.74303225 2.75866667
## std.dev coef.var
## 1.66092344 0.05548742
## ------------------------------------------------------------
## df[[1]]: SPF
## nbr.val nbr.null nbr.na min max range
## 9.00000000 0.00000000 0.00000000 30.90000000 38.30000000 7.40000000
## sum median mean SE.mean CI.mean.0.95 var
## 308.40000000 34.90000000 34.26666667 0.82276634 1.89730257 6.09250000
## std.dev coef.var
## 2.46829901 0.07203207
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.5065 0.4892
## 13
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.000141 0.0001413 0.167 0.689
## Residuals 13 0.011001 0.0008462
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 67.60 67.60 14.05 0.00243 **
## Residuals 13 62.53 4.81
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 4.198e-05 4.198e-05 9.805 0.00795 **
## Residuals 13 5.566e-05 4.280e-06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.000141 0.0001413 0.157 0.699
## df[[3]] 1 0.000210 0.0002098 0.233 0.638
## Residuals 12 0.010791 0.0008992
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.058262 1 64.7902 3.531e-06 ***
## df[[1]] 0.000349 1 0.3880 0.5450
## df[[3]] 0.000210 1 0.2334 0.6377
## Residuals 0.010791 12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
